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INTRODUCTION 


System identification is the process of constructing a mathematical model from input and output 
data for a system under testing, and characterizing the system uncertainties and measurement 
noises. The mathematical model structure can take various forms depending upon the intended 
use. The SYSTEM/OBSERVER/CONTROLLER IDENTIFICATION TOOLBOX (SOOT) is a 
collection of functions, written in MATLAB ' language and expressed in M-files, that implements 
a variety of modern system identification techniques. For an open-loop system, the centtal 
features of the SOCIT are functions for identification of a system model and its corresponding 
forward and backward observers directly from input and output data. The system and its 
observers are represented by a discrete model. The identified model and observers may be used 
for controller design of linear systems as well as identification of modal parameters such as 
dampings, frequencies, and mode shapes. For a closed-loop system, the central features of the 
SOCIT include identification of an open-loop model, an observer and its corresponding 
controller gain directly from input and output data. The basic package is capable of: 


1- Identifying system, forward and backward observer Markov parameters (pulse responses) 
from input and output time histories. 

2- Constructing a state space model from pulse responses. 

3- Identifying a state space model and its corresponding forward and backward observer gains 
directly from input and output time histories. 

4- Identifying a forward observer/Kalman filter gain with a given state space model, and input 
and output time histories. 

5- Computing variance and bias for identified modal parameters using Monte Carlo and 
perturbation procedures. 

6- Computing forward prediction errors and backward smoothing errors for any of the models 
generated. 

7- Identifying a state space model, and its corresponding controller gain and observer/Kalman 
filter gain directly from input, output and control force time histories. 

The unique features of this package are: 

1 - No nonlinear programming involved. 

2- No a priori noise information required. 

3- Guided model order selection. 

4- Direct identification of system & observer/Kalman filter. 

5- Direct identification of closed-loop controller. 

6- Suitable for stable & unstable systems. 


1 © Copyright 1985-91, by Mathworks, Inc. All rights reserved 
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SYSTEM/OBSERVER/CONTROLLER 
IDENTIFICATION TOOLBOX 

Reference 


Page 

arx b - calculates backward observer Markov parameters and residual error 1 1 

arx bat - calculates observer Markov parameters and prediction error 1 1 

arxc - computes the combined observer/controller Markov parameters 

from feedback control input and output data 14 

arx fb - calculates forward and backward observer Markov parameters, 

and their residual errors 11 

arx ps - calculates observer Markov parameters from pulse response samples 16 

bk_diag - transforms the modal form into a real block diagonal form 18 

block tr - computes matrix block transposition 20 

cpulse - converts rich input responses to pulse response time histories 21 

era - identifies a state space model from pulse response time histories 

using system realization theory (ERA) 23 

eradc - identifies a state space model from pulse response time histories 

using a data correlation technique (ERA/DC) 23 

freq_pt - plots the transfer function representation of a discrete time system 27 

hankl - forms a Hankel matrix from Markov parameters 29 

hankldc - forms a data correlation matrix for eradc 30 

k_abcd - identifies an Observer filter gain matrix from test data for a known 

discrete system model to whiten the stochastic residual. 31 

m2p - rearranges a Markov parameters sequence in the form of pulse 

response samples 33 

mar com - computes a specified number of system Markov parameters from 

observer Markov parameters 34 

mar_OC - computes a specified number of system, observer, and controller 

Markov parameters from observer/controller Markov parameters 36 

mar_sep - equivalent to mar_com but with separate observer Makov parameters 34 

mar yoc - computes a specified number of system, observer, and controller 

Markov parameters from feedback control inputs and outputs 38 

match - matches the eigenvalues identified from forward and backward models 40 

modal - computes a reduced stable or unstable model in modal coordinates 42 

monera - calculates variance of ERA identified parameters using the Monte 

Carlo approach 43 

OCid - identifies a state space model, an observer gain, and a controller 

gain from closed-loop experimental data 45 


3 



okid 

okidb 

okid_fb 

okidp 

p2m 

peradc 

prederr 

pred_efb 

prederb 

pulse 

ryucovar 

separate 

svpm 

svra 

uy_stack 

y_closed 

yesti 

y_pred 

yucovar 

yucovfb 

yucovb 

yycovar 


- identifies simultaneously a state space model and an observer gain 

from input and output data 50 

- modified version of okid using a backward observer 50 

- combined version of okid and okid_b using forward and backward 

observers 50 

- identifies a state space model from pulse response samples 50 

- rearranges pulse response samples in the form of Markov parameters 

sequence 62 

- calculates the variance and bias of ERA/DC identified parameters 

for single input and single output systems 63 

- computes prediction error from estimated observer Markov 

parameters 65 

- computes prediction and smoothing errors from estimated 

forward and backward observer Markov parameters 65 

- computes smoothing errors from estimated backward parameters 65 

- computes pulse response samples from general input and output data 67 

- computes the left correlation matrix associated with the feedback 

control input for an observer/controller identification 69 

- separates a given matrix sequence into two sequences 71 

- calcalates modal observability matrix and the singular value 

contribution of each mode to the pulse response samples 72 

- identifies a state space model from input and output data using a 

state vector realization technique 74 

- computes a stacked matrix with inputs and outputs 76 

- reconstructs closed-loop response time histories using ocid 

identified system, observer, and controller gain matrices 77 

- reconstructs outputs using an identified observer 78 

- reconstructs outputs using an identified system model 79 

- computes auto correlation and cross-correlation matrices between 

inputs and outputs 80 

- computes forward and backward auto-correlation and cross- 

correlation matrices between inputs and outputs 80 

- computes backward auto-correlation and cross-correlation matrices 

between inputs and outputs. 80 

- computes the left and right output residual correlation matrix 83 
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Road Map for Identification of an 
Open-Loop System 

(general input/output data) 


okid_b 
okid fb 



Note that Markov parameters also mean pulse response time histories. 
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Road Map for Identification of an 
Open-Loop System 

(pulse response time histories) 


okid_p 



Note that Markov parameters also mean pulse response time histories. 
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Road Map for Identification of a 
Closed-Loop System 

(general input/output data) 


ocid 



Note that Markov parameters also mean pulse response time histories. 
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Observer/Kalman Filter Identification (OKID) 

(Hubble Space Telescope Flight Test Data) 
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Closed-Loop Observer/Controller Identification (OCID) 

(Aircraft Flutter Test Data) 
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Purpose: 


arx_b,arx_bat,arx_fb 


Compute observer Markov parameters. 
Synopsis: 


[Ybf\=arx_b(m,r,u,yj),icl) 

[Yf\=aixJ)'dt(m,r,u,yj),ic2) 

[Yf,Yb]=aTx_fb(m,r,u,y t p,ic2) 

Description: 


The function computes observer Markov parameters from input/output data. The identified 
observer system is deadbeat of order p. Given l samples, r inputs, and m outputs, the input 
matrix u must have dimensions / x r and output y Ixm. Multiple experiments may be 
used in these functions. In that case, the input matrix u becomes / x ( rn, ) where n, is the 
number of experiments, and the output matrix becomes lx(mn t ). Function arx_bat solves 
the least squares problem of a forward observer; 


where 




y,=[y(0)y(l)-y(/-l)] 
Y ( = [D CB CAB — 


Vs* 


u( 0) m(1) u( 2) 
v(0) v(l) 

v(0) 


CA p - r B] 


u(l-D 

v(l- 2) 

v(/-p-l)J 


v(i) = U( ‘l ; * = 0,1 /-: 

LKOJ 


Here A = A + GC, B =[B + GD -G], in which A, B, C, D are system matrices, and G is 
the forward observer gain matrix. See function okid for more discussion on the definition 
of these matrices. The solution is stored in Yf of dimension mx[(r + m)p + r] If the 
experiment started from rest ic2(l)=0, otherwise, ic2(l)- 1. Once an estimate of the 
parameters is available the user is given the option to compute the prediction error 


e f -y f - Y f v f 


This computation when analyzing long records is time consuming. The square root of the 
diagonal elements of the inverse correlation matrix are proportional to the parameter 
variance. A chart depicting these values is plotted along with the prediction error. To bypass 
the prediction error option, set the second element of ic2 to one, i.e. ic2(2)=0. 

Function arx_b solves the least squares problem of a backward observer as follows 


where 


Z, = M 


= MO)y(l)— y(l-p-2)] 

Y b =[D + CB CA P ~'[GD — gJ CA p 2 [AB + GD ~g] 

r«(0) «(1) u( 2) - u{l-p-2) 

v(p) v(p + 1) v(p + 2) — ' v(/-2) 


v(l) v(2) v(3) v(l- p- 1) 


v(/) = 



i = 0,l. 


,/-l 


c[ab+gd -G]] 


Here, A = A~ x + GC, B = - A~ X B , in which /l, B, G, D are system matrices, and G is the 
backward observer matrix. See function okid_b for more discussion on the definition of 
these matrices. The solution is stored in Yb of dimension mx[(r + m)p + r]. Once an 
estimate of the parameters is available the user is given the option to compute the smoothing 
error, instead of the prediction error as in the case for the forward observer. 


The square root of the diagonal elements of the inverse correlation matrix are proportional to 
the parameter variance. A chart depicting these values is plotted along with the smoothing 
error. To bypass the smoothing error option, set the the variable ic to one, i.e. icl=0. Note 
that id is a scalar whereas ic2 is a vector with two elements. 

Observation of the forward and backward formulations shown above immediately reals that 
one may simultaneously compute forward and backward observer parameters. Both 
matrices are very similar in the sense that their lower sides are identical. Function arx_fb 
solves for backward and forward observer parameters simultaneously. All parameters used 
above also apply to this function. 

Algorithm: 

First, the correlation matrices are computed without actually constructing the individual 
matrices. The parameter estimate is obtained by 

Y = yV T (W T ) + 

where ( + ) refers to pseudo inverse. The pseudo inverse i$ computed using singular value 
decomposition. 
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Norm. Pred. Error 


Example: 


r=l;m=l;ic=[0 l|;index=0;p=2;L=100; 

a=[0 -0.16; 1-1]; b=|0 l]';c=[0 1]; d=0;G=[0.16 1]’; 

u=rand(L,m); 

y=dlsim(a,b,c,d,u); 

psize=r+p*(r+m); 

[Y]=arx_bat(m,r,u,y,p,ic); 

Compute Prediction Error (l=yes,0= no) =: 1 
Square Fitting Error Normalized 
1.9097e-29 




Q 1. ... | ...1 A J A .,-1 J A J 

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 

Parameter Number 


Y = [0.0000 1.0000 -1.0000 0.0000 -0.1600} 

See also: 

okid, okid_b, okid_fb 
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arxc 

Purpose: 

Compute combined observer/controller Markov parameters. 

Synopsis: 

[Ybar]=aT\c(y,ufb,ue,p, truncate) 

Description: 

This function computes combined observer/controller Markov parameters Ybar from 
feedback control input ufb, additive input excitation ue, and closed-loop response y. 

Consider a linear discrete system of the form 

x(k + \) = Ax(k) + Bu(k) 

y(k) = Cx(k) + Du(k) 

which is operating in closed-loop. The input to the system, u(k), consists of the 
feedback signal u^ijk) provided by an existing state feedback controller with gain F, 
and an additive excitation input signal u e (k) 

u(k) = u /b (k) + u,(k) 

- -Fx(k) + u t (k) 


The estimated state x(k) is provided by an existing observer of the form 

x (k + 1) = Ax(k) + Bu(k) - G d \y(k) - y(*)] 

y(k) = Cx(k) + Du(k) 

The output y(k) is the system closed-loop response due to an excitation u e (k). The 
function arxc solves for the observer/controller Markov parameters in Ybar which 
consists of 


D" 

0 


and 


;. c j 


(. A + GCy~'[B + GD -G] , k = 1, 2, .... p 


where A, B, C, D, and F are of the closed-loop system in operation, and G is another 
observer gain for the system such that 


C' 

-F 


( A + GC) k ~'[B + GD 


- G] = 0 , 


k = p + l, p + 2, ... 


where p is a number specified by the user. The number truncate specifies the number 
of data points to be deleted prior to application of the algorithm. This value is equal to 
the number of time steps that is expected for the existing observer to converge. 



Example: 

An example data file is contained in the file xsamp712. 
load xsamp712 

ue( 1 :600)=[l;y(l :600)=[] ;u( 1 :600)=[]; 
[Ybar]=arxc(y,-u,ue,20,300) 

Ybar = 


Columns 1 through 7 


-0.0932 0.2513 0.8236 
0.1228 -0.0551 -0.7010 


Columns 8 through 14 

0.0686 0.0149 -0.2200 
0.0247 -0.1287 0.0462 


Columns 15 through 21 

-0.2484 0.0971 0.0189 

0.0932 -0.0401 0.0931 

Columns 22 through 28 

0.1137 -0.0833 -0.0444 
-0.0387 0.0449 -0.0070 

Columns 29 through 35 

0.0894 0.1423 -0.0596 
-0.0062 -0.0320 -0,0272 

Columns 36 through 41 

0.0041 0.0171 -0.0598 

0.0579 -0.0791 0.0609 


-0.0771 0.3786 

0.0188 -0.1409 


0.2012 0.1215 
-0.0793 -0.0044 


-0.0101 0.0701 

-0.0382 0.0814 


-0.0189 0.0115 
0.0210 -0.0074 


-0.0859 -0.0226 
-0.0006 -0.0963 


0.0512 0.0595 

-0.0676 0.0661 


-0.0357 -0.5592 
0.0264 -0.0888 


-0.0940 -0.1043 
0.0364 -0.0138 


-0.1114 -0.0672 
-0.0203 0.0687 


-0.0695 -0.0910 
0.0147 0.0149 


0.0041 -0.0683 
0.0263 -0.0956 


-0.0819 

-0.0833 


Algorithm: 

The observer/controller Markov parameters are computed from feedback control input, 
additive excitation input, and closed-loop response data. These parameters are used in 
the function ocid to compute a realization of the system state space matrices, the existing 
controller gain, and an observer gain. 

See also: 

ocid, mar_yoc, mar_oc, ryucovar, y_closed, separate 


arxps 

Purpose: 

Compute observer Markov parameters from pulse response histories. 

Synopsis: 

\d,ys ,yo]=aTx_ps(y,m,p,ic) 

Description: 


The function computes observer Markov parameters from pulse response time histories. The 
identified observer system is deadbeat of order p. Given / samples, r inputs, and m outputs, 
y is / x mr .The least squares problem solves the following equations 

y = YV 

where 


y = bi y 2 — y r l >,0) — y,(/-i)];/ = i,...,/-. 

f = [K, y 2 ] ; y;=[D c(b+gd) ca(b+gd) ... ca p ~\b+gd)i 

Y 2 =[-CG -CAG - -CA p -'G]; 


v = [v, v 2 



II 

\(0) u,( 1) 
«,(0) 

«,( 2) ••• 
U,(l) ... 

r _ 

• i 
=sf "5* 

K» 

II 

w S' 
© 

0 

1 

- y,(/-2)’ 

- x(/-3) 



Mj(0) ... 



L : 



0 


M,(0) = 


1 ( ith element ) 


0 


;u,(*) = 0; k = 


The solution is stored in [d, ys, yd\ where d is the system transmission matrix D, 

ys = [C(B + GD) CA(B + GD) ... CA P '\B + GD)] 
and 

yo = [-CG -CAG - -CA P ~'G] 

The matrix ys has dimension mxrp and yo has dimension m'xmp. 

Algorithm: 

First, the correlation matrices are computed without actually constructing the individual 
matrices. The parameter estimate is obtained by 
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Y = yV T (W T ) + 

where ( + ) refers to pseudo inverse. The square matrix W T has the following special form 


VV T = 


I a 

1 r(p+])xr(p+\) U'rip+Vxmp 


a 


a 


mpx(p+\)r Pmpxmp J i=l 

0.„ b,<0) - ?,«))] [>,(1) - *,(!)] 

[>,(0) - >,(0)] 


wv>x(p+l)r 


bi(p-i) — y,(p-i)] 
bi(p-2) — y r (p- 2)] 

bx (0) - y,(0 )] 


which make the pseudo inverse (W r ) + easier as shown below 


(W T ) + 


I + a T ((3-aa T y a -cc T (f}-aa T y 
~(P~aa T ya (fi-aa 7 )* 


The pseudo inverse (/? - aa T Y is computed using singular value decomposition. Once an 
estimate of the parameters is available the user is given the option to compute the prediction 
error 

e = y-YV 

This computation when analyzing long records is time consuming. The square root of the 
diagonal elements of the inverse correlation matrix are proportional to the parameter 
variance. A chart depicting these values is plotted along with the prediction error. To bypass 
the prediction error option, set ic=0. 

Example: 

This example is to identify observer Markov parameters from the pulse response samples of 
a three-mass- spring-dashpot system with two inputs and one output. 


k 1 = 1 .0;k2=2.0;k3=3.0;ml = 1 .0;m2= 1 ,0;m3= 1 .0;ratio=2*0.005; 

K=fkl+k2 -k2 0 ; -k2 k2+k3 -k3; 0 -k3 k3]; 
Khalf=sqrtm(K);Damp=ratio*Khalf; 

Ac=[zeros(3,3) eye(3,3);-K -Damp];Bc=[zeros(3,2);l 0; 0 1; 00];C=[zeros(l,5) 1]; 

dt=1.0;pt=100;p=10;[m,n]=size(C);[n,r]=size(Bc);D=zeros(m,r); 

t=[dt:dt:pt*dt]'; 

[A,B]=c2d(Ac,Bc,dt); y=[]; 
for i=l:r; 

y=[y dimpulse(A,B,C,D,i,pt)l; 
end; 

[d,ys,yo]=arx_ps(y,m,p,0) 

H=mar_sep(ys,yo,d,m,r, 1 0) 

See also: 

mar_sep, okid_p, arx_bat, arx_ps 


1 7 


bkdiag 


Purpose: 

Transform the complex modal form of a discrete model into a real block diagonal form. 
Synopsis: 

[Ab,Bb,Cb]=bk_d\ag{Lambda,B,C) 

Description: 

Given a complex diagonal model 


x(k + l) = Ax(k) + Bu(k); 
y(k) = Cx(k ) + Du(k) 


where 

A = diag(X j,..., X t , ct s+ 1 + ■ ffj+i — jfii+i + jf$ H , a H — jfl n ) 

' Ih 

b, 

B= a s+l + JP,+l 

ttj+l ~ jflt* 1 
+ jPn 

, <**~jP* 

which is a complex model, the function bk_diag transforms this model to the following 
block-diagonal form 

. x h (k + 1) = A b x b (k) + B b u(k)\ 

y(k) = C b x b (k) + Du(k) 

where 


;C = [c, ••• c s r l, +l +jn, +l r ?, +1 


X 




«, + . & + i 

-Ps+I a, + i 

«. P„ 
"A. 
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A 


B b = 


2a 

-2A 


s + 1 


♦ Q — [ f 


A 


*♦1 


‘*+1 


/*„] 


2a, 

-2ft _ 


All the variables in this block-diagonal form are real rather than complex as in the diagonal 
form. This function is used in conjunction with function modal to reduce an identified 
model (stable or unstable) to a stable real block diagonal model for numerical simulations to 
compare with real data. 


Example: 

randCnormal'); 

n=7; 

am=rand(n,n);bm=rand(n,2);cm=rand(l,n); 
[v,lambda]=eig(am); 
lambda=diag(lambda); 
bm=vM>m; cm=cm*v; 
flambda,k]=sort(lambda); 
bm=bm(k,:);cm=cm(:,k); 

[ a,b,c 1 =bk_diag(lambda,bm,cm) 


1.0188 0 0 
0 -0.6935 -1.1752 

0 1.1752 -0.6935 

0 0 0 
0 0 0 
0 0 0 
0 0 0 
b = 

-3.5006 -0.9259 
-0.2835 1.3741 

1.6978 2.7562 

3.8881 0.6494 

-1.5251 -2.3995 
-3.3324 -0.0417 
2.1830 0.0425 


0 

0 

0 

0.4642 

1.8936 

0 

0 


0 0 0 

0 0 0 

0 0 0 

-1.8936 0 0 

0.4642 0 0 

0 1.9055 -0.8629 

0 0.8629 1.9055 


C -0.0730 -0.0421 -0.2797 -0.2737 0.3682 -1.0026 -0.4538 
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See also: 

okid_b. okid_fb, modal 


b!ock_tr 

Purpose: 

Compute the matrix block transposition. 

Synopsis: 

[at]=b\ock_tir(nrow,nblock,ncol,aJlag) 

Description: 

The function performs a 2-way block transposition depending on the input matrix. When 
flag is set to 1 the wide matrix 


is block transposed to 


a = [>’o A -W*] 




L-^ nblock J 

where each block yj is a matrix of dimension nrow x ncol. The reverse operation is obtained 
when inputing a tall matrix and flag = 0. 


Example: 

y0 = [0 1; 2 3]; 
yl =[4 5; 6 7|; 
a=(y0yll 

a- 

0 14 5 

2 3 6 7 

[at|=block_tr(2,2,2,a, 1 ) 

at = 

0 1 
2 3 

4 5 

6 7 

a=block_tr(2,2,2,at,0) 

a- 

0 14 5 

2 3 6 7 

See also: 


mar_com 
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cpulse 


Purpose: 

Compute pulse response samples using FFT. 

Synopsis: 

(y.v]=cpulse(y,u,/',m); 

Description: 

The function cpulse calculates the unit pulse response samples (Markov parameters) from 
input and output time histories by using a frequency domain approach. The input time 
histories are required to be sufficiently rich (e.g. random inputs). The system input and 
output histories u(t) and y(t) must be stored as follows 


u = 


y = 


m„( 0) 
u uO) 
Mil (2) 

M,,(/-l) 

yn(O) 

y.i0) 

>-1.(2) 


u r \ ( 0 ) 
M r ,(D 

«rl(2) 

u rl (l~ 1) 

y-iO) 

y«i(2) 


M,j(0) 

M U (1) 

Mu(2) 

«u(/-l) 

y i.0) 
y i.(2) 


M„(0) ‘ 

M rj (l) 

M rj (2) 

M„(/-l) 

y~«>) 

>wd) 

y«(?) 


y u (/-i) ••• y m i(/-i) y\ jd _1 ) ^ 


where r is the number of inputs and m is the number of outputs, and M >> (t)(y i/ d)) « the i-th 
input (output) of the )-th test at discrete time t. The number of experiments s should be 
greater or equal to r, integer / is required to be even. The ouput of this function is the pulse 
response histories ys. 


Example: 

The example is to compute the pulse response samples from 5 data sets of random inputs for 
a single-input and single-output second-order system with sampling interval dt=0.2, pt-MZ 
samples, natural frequency co„ = 1 and damping factor £-0.1, 

a=[0 l;-l -0.2];b=l0; l|;c=[l 0];d=0; 

pt=5 1 2;dt=0.2;rand(’normal'); 

t=[dt:dt:pt*dt]'; 

u=rand(pt,5); 

for i=l:5, 

yt(:,i)=lsim(a,b,c,d,u(:,i),t);end 
[yl]=cpulse(yt,u,l,l); 
plot(yl),title('Pulse response') 
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Algorithm: 


The function cpulse uses a frequency domain approach to compute the pulse response 
samples (Markov parameters) from output time histories generated from rich input signals. 
First the FFT is applied to calculate the discrete frequency response functions of input u(t) 
and output y(t). Then the input and output frequency response functions, are used to 
compute the discrete transfer functions G(z), and the pulse response samples are the inverse 
FFT of G(z). 

Y(z)=G(z)U(z), Y(t)=FFT-'(G(z)) 


See also: 
pulse 
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era,eradc 

Purpose: 

Identify a state-space model from pulse response samples (Markov parameters). 

Synopsis: 

[a,b,c,d,sg,eg,mh]=era(y,m,r,n,nm,dty, 

\a,b,c,d,sg,eg,mh]=cT&dc(y,m,r,n,nm,dt); 

Description: 

The function era identifies a state-space model of a multi-input and multi-output linear, time- 
invariant system from pulse response samples (Markov parameters). The pulse response 
samples arc stored as 


' >11 ( 0 ) 

>ml(°) 

>lr(0) 

>w(0) 

>11 ( 1 ) 

- >ml 0) 

>lr0) 

— >mr( 1) 

_>,l(/-l) 

- >m 

... > lr (/ — 1) 

- > W (/- 1) 


where y tj {t) is the i-th output at discrete time t due to a unit pulse at the y'-th input. The 
system to be identified has r inputs and m outputs. The identified model order is chosen as 
n. When n is set to zero, user's inputs is required on-line to specify the desired order of an 
identified model, based on the singular values of a Hankel matrix. Scalar nm specifies the 
number of sample shifts for forming the rows of a Hankel matrix (see the Algorithm section 
few definition). Integer nmxm should be greater than the model order n. Scalar dt specifies 
the data sampling interval. [ A,B,C,D,sg,eg,mh]=CT&(y,m,r,n,nm,dt ) returns an n-th order 
linear, time-invariant identified discrete model: 

x(k + l) = Ax(k) + Bu(k) 

y(k) = Cx(k) + Duik) 

Matrix eg contains modal parameters of the identified model with damping ratios (%) in the 
first column and frequencies (Hz) in the second column. The third column of eg gives the 
eigenvalues of the corresponding continuous-time model. Note that the identified discrete 
model can be easily transformed to a continuous-time model. Vector sg , whose elements are 
singular values of the Hankel matrix, can be used as reference to choose the model order n. 
The first column of matrix mh gives the normalized singular contribution of each identified 
mode in matrix eg to the pulse response samples whereas the second column gives the 
modal amplitude coherence. The maximum singular value is chosen to normalize the first 
column of matrix mh. These normalized singular values are used to weight the importance 
of the individual mode to the pulse response samples. Each element in the second column 
of mhi (the i-th element of mh) is between 0 and 1; mh t — > 1 indicates that the identified 
mode eg; (the i-th eigenvalue) is reliable. 

The function eradc is similar to era, but it uses the ERA data correlation method [Juang88] 
to identify the system model. For a long data lenght, eradc is recommended for use because 
it takes less memory and computational time to solve for singular values of the Hankel 
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matrix. In eradc, the size of the data correlation matrix has been minimized to save time in 
computing its singular values. 

Examples: 

Example 1 : 


i) Calculate pulse response samples of a single-input and single-output second-order 
system with sampling interval dt=0.2, natural frequency G>„ = 1, and damping factor 
6 = 0 . 1 . 

ii) Use era interactively to identify a model from the pulse response samples from (i) and 
transform the identified discrete model to a continuous-time model. 

iii) Plot the error between the pulse response samples from (i) and the pulse response 
samples of the identified model. 

iv) Display eigenvalue matrix eg, singular value vector sg and modal amplitude coherence 
matrix mh. 

a=[0 1;-1 -0.2]; b=[0;l];c=[l 0];d=0; 

pt=100;dt=0.2;t=[dt:dt:pt*dt]'; 

u=zeros(pt,l);u(l,l)=1.0; 

y=lsim(a,b,c,d,u,t); 

clg 

[a 1 ,b 1 ,c 1 ,d 1 ,sg,eg,mh]=era(y, 1 , 1 ,0,20,dt); 

y 1 =dlsim(a 1 ,b 1 ,c 1 ,d 1 ,u); 

eiror=y-yl; 

clg 

plot(error),title('Pulse response error') 
pause 

eg,mh,sg(l:10) 

Example 2: 

i) Calculate pulse response samples from 6 data sets for a single-input and single-output 
second order system with sampling interval dt=0.4, /=512 samples, natural frequency 
0) K = 1, damping factor £ = 0.1 , and noise standard deviation 0.04. 

ii) Use eradc interactively to identify a model from the pulse response samples from (i) 
and transform the identified discrete model to a continuous-time model. 


iii) Display eigenvalue matrix eg, singular value vector sg, and modal amplitude coherence 
matrix mh of the eradc identified model. 

iv) Plot the pulse response samples of the original model and the eradc-identified model. 
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a=[0 1 1 -0.2]; b=[0;l];c=[l 0];d=0; 

dt=0.4;pt=5 1 2;t=[dt:dt:pt*dt]’; 

randCnormal'); 

u=rand(pt,6); 

for i=l:6, 

yt(:,i)=lsim(a,b,c,d,u(:,i),t);end 

y=yt+0.04*rand(pt,6); 

yi=cpulse(y,u,l,l); 

clg 

[a 1 ,b 1 ,c 1 ,d 1 , sg,eg,mh] =eradc(y i, 1 , 1 ,0,50,dt); 

eg,mh,sg(l:10) 

clear u 

u=zeros(pt, 1 );u( 1 , 1 )= 1 .0; 

y 1 =dlsim(al ,bl ,c 1 ,d 1 ,u); 

yO=lsim(a,b,c,d,u,t); 

y=ry0yl]; 

clg 

plot(y),title(’Pulse response') 
pause 


Algorithm: 

The function era uses the Eigensystem Realization Algorithm (ERA) from [Juang85], which 
uses Markov parameters (pulse responses) to form the Hankel matrix 



' Y J 

Y j + 1 

1 

+ 

H(j-\) = 

Y j+t 

V 

■” Y j+P+ 1 
• « 

m « 

« • 


/j+r 

^/+ y+ 1 

"■ Y j+y + p_ 


>ll(0 •" y\r(i) 

^2l(0 •” Ji2r(0 


bml(‘) - ymr 0')J 

where Y t is the i-th Markov parameter and y (> is the i-th output at discrete time t due to a unit 
pulse at the j - th input. From the measurement Hankel matrix, ERA uses the SVD of H(0), 
H( 0) = U X V T , to identify a Jfc-th order discrete state-space model as 

At = r t m ulH(\)V k -Li' n 
B t = li ,2 V, r £, 

Q = E T m U k £i' 2 
D t = Y( 0) 

where matrix X* is the upper left hand k x k partition of X containing the k largest singular 
values along the diagonal. Matrices U k and V k are obtained from U and V by retaining only 
the k columns of singular vectors associated with the k singular values. Matrix E m is a 
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matrix of appropriate dimension having m columns, all zero except that the top m x m 
partition is an identity matrix. E r is defined analogously. 

The function eradc uses a special case of the ERA/DC algorithm from [Juang88]. It starts 
with the Hankel matrices H( 0) and 7/(1) to generate the block correlation matrices 
R(i) = 0) r . The SVD of /?(0) = t/£ V T , is applied to identify a kth order model as 

A t =ii U 2 ulKmv t ii m 
S*=X*' 2 V/W(0)E, 
c t = E T m U t si ' 2 
D t = Yi 0 ) 

The function eradc uses the block correlation matrices R(i)= H(i) T H( 0) if this matrix is 
smaller than /?(/) = H(i)H( 0) T . 

Limitations: 

The most time consuming step in each algorithm is the singular value decomposition. The 
number of floating point operations for SVD are roughly a cubic function of the matrix 
dimension. Also the SVD of a large matrix needs a lot of memory. The size of the SVD 
matrix in era and eradc is (nmxm)x(£- nm-lxr) and (nmxm)x(nmxm) [or 
( t-nm)x ( t-nm ) if this is smaller] respectively where i is the length of the data. In era, 
the column number of the Hankel matrix may be very large if the data length is large. 
However, the column number in eradc can be chosen as large as desired, because the size of 
the data correlation matrix depends only on the number of rows. 

See also: 

cpulse, okid, okid_b, okid_fb 
References: 

[1] Juang, J. N. and Pappa, R. S., "An Eigensystem Realization Algorithm for Modal 
Parameter Identification and Model Reduction," Journal of Guidance, Control, and 
Dynamics, Vol. 8, No. 5, 1985, pp. 620-627. 

[2] Juang, J. N., Cooper, J. E., and Wright J. R., "An Eigensystem Realization Algorithm 
Using Data Correlation (ERA/DC) for Modal Parameter Identification," Control Theory and 
Advanced Technology, Vol. 4, No. 1, pp. 5-14, 1988. 

[3] Lew, J. S., Juang, J. N. and Longman, R. W., "Comparison of Several System 
Identification Methods for Flexible Structures," Proceedings of the 32nd Structures, 
Structural Dynamics, and Materials Conference, Baltimore, MD, April 1991, pp. 2304- 
2318. 

[4] Juang, J. N., "Mathematical Correlation of Modal Parameter Identification Methods Via 
System Realization Theory," International Journal of Analytical and Experimental Modal 
Analysis, Vol. 2, No. 1, Jan. 1987, pp. 1-1 8. 
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freqpt 

Purpose: 

Plot the transfer function representation of a discrete time system. 

Synopsis: 

[mag phase} =frcq_pt(a, b,c,d,p,dt,iu) 

Description: 

Given a discrete model 


x(k + \) = Ax(k) + Bu(k) 

y(k) = Cx(k) + Du(k) 

the function computes the transfer function representation given by 

G(e JaAT ) = C{e ieAT l - AT' B + D 

The parameter p determines the number of spectral lines plotted, dt is the sample time, iu is 
the input number to be plotted. The function [mag,phase]=freq_pt(/4, B, C,D,p,dt,iu) 
returns the magnitude (mag) and phase plots of the transfer functions for the m-th input. 
The function prints the value for the Nyquist frequency. The user is prompted for lowest 
frequency to plot, and for the upper frequency bound. The upper frequency bound must be 
given in terms of percentage of the Nyquist frequency. These values are used to define the 
frequency range. The actual plot scale may be slightly different because it is determined by 
the plotting function. 


Example: 


a=[0 1;-100 -0.0021; b=[0;ll;c=[l 0];d=0; 

dt=0.04;pt=1024;t=fdt:dt:pt*dtj'; 

rand('normar); 

u=rand(pt,l); 

y=lsim(a,b,c,d,u,t); 

y=y+0.0042*rand(pt,l); %3 percent noise 
[a,b,c,d,m]=okid( 1 , 1 ,dt,u,y, 'batch', 10); 
[mag,phase] =freq_pt(a,b,c,d,300,dt, 1 ); 

Nyquist frequency (Hz) is =: 12.5 

Enter lower frequency to plot (Hz)=: 0.1 
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See also: 





hankl 

Purpose: 

Form a Hankel matrix from a sequence of Markov parameters. 

Synopsis: 

l£>,//]=hankl(M arkov.p) 

Description: 

Given a discrete model 


x(k + 1) = Ax(k) + Bu(k) 
y(k) = Cx(k) + Du(k) 

with r inputs and m outputs, the pulse response samples are typically stored as 

Y P =[y t ^ - xj 

where y t = response samples due to a unit pulse at the /- th input. Each y ( has m columns 
and l rows where t is the length of data. The pulse response samples can be rearranged 
by using function p2m to the following sequence of system Markov parameters 

n=[n y, n - n-0 

=[d cb CAB — CA'-’b] 

Function [D,//]=hankl(K5',p) return the transmission matrix D and a Hankel matrix defined 


Y <- 2 


H( 0 ) = 


Y. Y. 


p+\ 


Note that all the data pass into the function are used to form the Hankel matrix. The inner 
products are used in matrix multiplication to reduce computational time. The size of the 
matrix is px(l- p- 1). The longer the data length is, the larger the number of rows of the 
matrix becomes. This function is used in era for identification of system matrices. 

See also: 

era, eradc, m2p 
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hankldc 

Purpose: 

Form a data correlation matrix from a sequence of Markov parameters. 

Synopsis: 

[D,R]=hankl(Markov,p) 

Description: 

Given a discrete model * •' • 


x(k + 1) = Ax(k) + Bu(k) 
y(k) = Cx(k) + Du(k) 

with r inputs and m outputs, the pulse response samples are typically stored as 

Y p=hx y 2 % - y m ] 

where y,= response samples due to a unit pulse at the i-th input. Each y l has m columns 
and t rows where t is the length of data. The pulse response samples can be rearranged 
by using function p2m to the following sequence of system Markov parameters 

n = [n y, n - 

= [D CB CAB - CA l ~ 2 B ] 

From this sequence of Markov parameters, define two Hankel matrix as 



m Y t n - Y t -„-2 



r Y Y ••• Y 1 

1 2 1 t~p-2 






K K ••• K 


1! 

i 

» 

sc* 

... 

; // = 


*2 1 3 t-p-\ 






Y y ••• Ym * 



Y Y . Y t , 



1 P 1 p + 1 V-3 



l p p+l t - 3 



Y Y ... Y 

l p + 1 V+2 * t-2 



Function [D,/?]=hankldc(K.s,p) returns the transmission matrix D and a data correlation 
matrix R defined as 

R = HHj 

Note that all the data pass into the function are used to form the data correlation matrix. The 
size of this matrix is (p + l)xp which is independent of the length of the data. This 
function is used in eradc for identification of system matrices. 

See also: 

era, eradc, m2p 
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k abed 


Purpose: 

Identify an Observer/Kalman filter gain matrix from test data for a discrete model to whiten 
the stochastic residual. 

Synopsis: 

[G, Gmar]=k_&bcd(A ,B ,C ,D ,u,y ,p) 

Description: 

The function k_abcd solves for an observer/Kalman filter gain matrix of a system in the 
form 


x(k + 1) = Ax(k) + Bu(k) - G\y(k) - y(*)] 
y(k) = Cx(k) + Du(k) 


where x(k) is the estimate of the state x(k) and y(k) is the estimate of y(k). The system has 
m outputs, r-inputs, and time samples dt apart. The input/output time histories are stored as 
column matrices. Initially an estimate of the desired observer Markov parameter number p 
must be given for whitening the residuals. It is suggested that p is chosen such that the 
product pxm is greater than or equal to the order of the system. The identified matrix, G, 
and observer gain Markov parameters, Gmar (stored as a column matrix), are returned to the 
main program. See references for detailed information. To identify a stable observer which 
whitens the residual between the real output y(k) and the estimated output y(k), the system 
matrices A, B,C, D, must be reasonably close to the true ones. 


Example: 


From the test data of a truss structure, use the function to compute a set of system matrices 
and an observer gain matrix. The order of the system is determined automatically by the 
function eradc. The residual is further whitened by using the function k_abcd to modify the 
observer gain matrix. The function y_esti is used to compute the estimated outputs which 
are then subtracted from the real output to obtain the residual. The following is output taken 
from a typical run. 


load sample 
fpt,i]=size(u); 

fa,b,c,d,ml=okid(n,r,dt,u,y,’batch',20); 
ftime,y_e]=y_esti(a,b,c,d,m,u,y,dt,pt); 
[ml,Cake]=k_abcd(a,b,c,d,u,y,20); 
[time,y_el ]=y_esti(a,b,c,d,ml ,u,y,dt,pt); 
clear a b c d u 


rey=y-y_e; 

[yvt, vvtl=yycovar(n,2,pt,rey, 1 ); 
clear res 


wt 
wt = 

2.1068e+02 

-2.5208e+01 

1.3304e+02 

7.4350e+00 


-2.5208e+01 

1.6014e+02 

1.6889e+01 

7.4626e+01 


1.3304e+02 

1.6889e+01 

2.1061e+02 

-2.5116e-K)l 


7.4350e+00 
7.4626e+01 
-2.51 16e+01 
1.6015e+02 


res/=y-y_el; 

[yvt, vvtl]=yycovar(res/,2,l); 


clear res 1 
wtl 
vvtl= 

1.2454e+02 -4.8890e+fil 
-4.8890e+01 1.3218e+02 
-1.9179e+00 2.6610e-01 
-2.999 le+OO -1.4728e-K)l 


-1.9179e+00 -2.999 le+OO 
2.6610e-01 -1.4728e+01 
1.2469e+02 -4.8836e+01 
-4.8836e+01 1.3220e+02 


The matrix 


wt = 


T res,(l)' 


res 2 ( 1) 
res, (2) 
res 2 ( 2) 

[res,(l) res 2 ( 1) res, (2) res 2 (2)] 


is the expected value of the auto- and cross-correlation of the residual obtained from the 
okid identified observer, whereas 



>esl,(l)" 


wtl = E 

resl 2 (l) 
resl, (2) 

[resl,(l) resl 2 (l) resl, (2) resl 2 (2)] 


resl 2 (2) 



is obtained from the k_abcd identified observer gain using the okid identified system matrices. 
It is obvious that the residual resl is whiter than res. 

Algorithm: 

The algorithm used here is similar to that for the okid which identifies a set of system 
matrices. A, B, C, D, as well as an observer gain matrix, G. For given system matrices A, 

B, C, D, the deterministic component of the output can be subtracted out. The observer 
gain is obtained by whitening the remaining residual. For details see references. 

§ee also: 

okid, okid_b, okid_fb, yycovar 

References: 

[1] Chen, C. W., Huang, J.-K., and Juang, J.-N., "Identification of Linear Stochastic Systems 
Through Projection Filters," presented at the AIAA 33rd Structures, Structural Dynamics & 
Materials Conference , 1992, Paper No. AIAA-92-2520. 

[2] Juang, J.-N., Chen, C. W., and Phan, M., "Estimation of Kalman Filter Gain from Output 
Residuals" NASA Technical Memorandum TM- 1 07603, Langley Research Center, 
Hampton, VA., March 1992. 
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m2p 

Purpose: 

Rearrange a Markov parameters sequence in the form of pulse response time histories. 
Synopsis: 

Y p =m2p(Y s ,r) 

Description: 

Given a discrete model 

x(k + \)=Ax(k)+Bu(k) 

y(k ) = Cx(k) + Du(k) 

with r inputs and m outputs, the sequence of system Markov parameters is defined as 

y,=[d cb cab ca 1 ~ 2 b] 

The pulse response samples are typically stored as 

Y p =[y i >2 - 3v] 


where y, = response samples due to a unit pulse at the i-th input. Each y i has m columns 
and / rows where t is the length of data. The sequences Y s and Yp are equivalent in the 
sense that both represent pulse response samples. The function converts Y s to Y p . Note 
that Y p is the sequence used in the eradc and era functions. 

If the combined system/observer Markov parameters (see function mar_com) are used, the 
input matrix B becomes \B G ] and the sequence Y s becomes 

Y S =[D C[B G ] CA[B G] — CA l2 [B G]] 

where the input number changes from r to r+m. The pulse response samples in Y p 
becomes 

y p = b, Yr y r+i y r +2 — ^+*1 

with m additional txm matrices due to the observer gain matrix G which is treated as an 
input matrix. The first r matrices represnt the system pulse response samples and the last 
m matrices mean the observer gain pulse response samples. 

See also: 

mar_com, okid, okid_b, okid_fb, p2m 


33 



Purpose: 


mar com, mar sep 


Recover the system Markov parameters from a set of observer Markov parameters. 
Synopsis: 

H=maT_com(ybar,r,njnarkov) 

H=mar_sep(yOb,y]b,T),r,n_niarkov) 

Description: 

Given an observer of the form 

*<* + l) = Ax(*) + [B 

y(k) = Cx(k) + Du(k) 

with r inputs and m outputs, the sequence of observer Markov parameters previously 
computed is passed to the function mar_com with 

[d c{b - G } ca{b -C} - C A P {B -G}] 

or to the function mar_sep with 

yOb = [CB CAB CA 2 B - CT^'b] 
and 

y\b = [CG CAG CliG ••• 

There is usually a small number, p, of nonzero observer parameters which is less than 
n^markov ,_7he observer matrices are related to the system matrices by 
A = A + GC, B = B + GD, and D is the direct transmission term which is the same for 
system and observer. The function computes recursively n markov parameters of the 
original system and puts them in the form 

H=[{D 0} C{B -G} CA{B -G} ••• CA P {B -G}] 

Algorithm: 

The parameters are computed using the following recursive formula 

CA k [B -G] = CA l B- Y i C~A‘GCA k , l [B -G]-CA*G[D 0] 

i=0 

for k=\,2,...n_markov 
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Example: 


a=[0 -0.16; l-l];b=[0 1|’; 
c=[0 1]; d=0; G=[0.16 1]'; 
abar=a+G*c; 
bbar=[b+G*d -G]; 
ybar=[d c*bbar c*abar*bbar]; 

[Hl=mar_com(ybar, 1 ,2) 
bg=[b -Gl; 

Hs=[d 0 c*bg c*a*bg] 

H = 

0 0 1.0000 -1.0000 -1.0000 0.8400 


Hs = 

0 0 1.0000 -1.0000 -1.0000 0.8400 

See also: 

okid, okid_b, okid_fb, okid_p 
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maroc 

Purpose: 

Compute a specified number of the system, observer, and controller Markov parameters 
from observer/controller Markov parameters. 

Synopsis: 

[H]=mar_oc(Ybar,r,m,Ntotal) 

Description: 

From a sequence of observer/controller Markov parameters arranged in Ybar in the 
following order 


D 

0 


, and 


C 

-F 


\(A + GC) k ~'[B + GD 


-G], * = 1,2, 


P 


the function computes the following sequence of system, observer, and controller Markov 
parameters 


D 

0 



~G], 


' C ' 
-F 


A[B — G] , ... 



~G] 


that are arranged in the matrix H in this order. The scalar r denotes the number of inputs, 
and m the number of outputs. Ntotal specifies the number of Markov parameters in H to 
be returned to the user. For background information, the sequence of observer/controller 
Markov parameters are computed from closed-loop excitation data, and the function mar_oc 
is then used to unscramble this sequence to obtain the system, observer, and controller 
Markov parameters. 

Example: 

load xsamp712 

ue(l:600)=[];y(l :600)=[ l;u( 1 :600)=[]; 

[Ybar]=arxc(y,-u,ue,30,300); 

[H]=mar_oc(Ybar, 1 , 1 ,30) 

H = 


Columns 1 through 7 


-0.0689 

0.1778 0.8012 

0.0572 

1.0144 

0.1181 

0.5774 

0.0351 

0.0422 -0.6740 

-0.1063 

-0.6234 

-0.0419 

-0.8527 

Columns 8 through 14 





0.0915 

0.4184 -0.1442 

0.2114 

0.0086 

0.0914 

-0.1890 

-0.1076 

-0.6742 -0.0692 

-0.6106 

0.0527 

-0.3847 

-0.0236 
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Columns 15 through 21 

-0.1854 -0.0041 -0.3831 -0.1030 -0.5743 -0.1404 -0.6633 

-0.1319 0.1225 0.2654 0.0310 0.5776 0.1159 0.8443 

Columns 22 through 28 

-0.0318 -0.7223 -0.0736 -0.7292 0.0206 -0.7679 -0.1022 

0.1240 0.9964 0.0591 1.0904 0.0730 1.1074 0.0075 


Columns 29 through 35 

-0.5967 0.0881 -0.4962 -0.0521 -0.2350 0.0872 -0.1258 
1.0717 0.0602 0.8314 -0.0597 0.5758 0.0144 0.2073 


Columns 36 through 42 

0.0442 0.0852 0.0671 0.2845 0.1206 0.5419 0.0896 
-0.1010 -0.0944 -0.0799 -0.4433 -0.1145 -0.7440 -0.1510 


Columns 43 through 49 

0.7482 0.1712 0.8393 0.0375 0.9105 0.0575 0.8713 
-1.0425 -0.1319 -1.2376 -0.1653 -1.2999 -0.0688 -1.2899 


Columns 50 through 56 

0.0530 0.8104 0.0221 0.6379 0.0741 0.3764 -0.0248 
-0.0444 -1.1347 -0.0149 -0.8797 0.0075 -0.5058 -0.0365 


Columns 57 through 61 

0.1169 -0.0111 -0.1443 -0.1061 -0.3409 
-0.0961 0.0786 0.3044 0.0666 0.7439 


Algorithm: 

The system, observer, and controller Markov parameters are computed from the 
observer/controller Markov parameters by a set of recursive equations. If more than p 
number of system, observer, and controller Markov parameters are required to be solved, 
the extra observer/controller Markov parameters are set to zero. For further details, see 
references. 

See also: 


arxc, mar_yoc, odd, ryucovar, separate 

Reference: 

[1 ) Juang, J. N. and Phan, M., "Identification of System, Observer, and Controller from 
Closed-Loop Experimental Data," Presented at the AIAA Guidance, Navigation, and 
Control Conference, Hilton Head, South Carolina, Aug. 10-12, 1992. 
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Purpose: 


Compute system, observer, and controller Markov parameters directly from feedback 
control input ufb, additive input excitation ue, and output response y. 

Synopsis: 

[ocs]=mar_yoc(y,ufb,ue,p,truncate,NtotaO 

Description: 

The function mar_yoc solves for the Markov parameters 


D~ 

0 



-G], 



A[B -G] , ... 



-G] 


from feedback control input ufb, additive input excitation ue, and output response y. The 
data is stored as column matrices. For example, for a system with m outputs, the closed- 
loop output data matrix y contains m columns and as many rows as the number of data 
points available. The data matrix ufb contains the feedback control signal, ue contains the 
additive excitation input signal. The number p denotes the number of observer/controller 
Markov parameters to be solved. The number truncate specifies the number of data points 
to be deleted prior to application of the algorithm. This value is equal to the number of time 
steps that is expected for the existing observer to converge. Ntotal is the total number of 
Markov parameters to be solved for from p identified observ er/c ontroller Markov 
parameters. This is done by setting the extra observer/controller Markov parameters to be 
zero. 



-G] = 0 , k = p + \, p + 2, ... 


These Markov parameters are in the form ready to be used to obtain a state space realization of 
the system matrices, observer gain, and controller gain. 

Example: 

An example data file is contained in the file xsamp712. 
load xsamp712 

[ocs]=mar_yoc(y,-u,ue,3,200,4) 


ocs = 


Columns 1 through 7 


-0.0689 

0.0351 


0.1778 

0.0422 


0.8012 

-0.6740 


0.0572 

-0.1063 


1.0144 

-0.6234 


0.1181 

-0.0419 


0.5774 

-0.8527 


Columns 8 through 14 


0.0915 0.4184 -0.1442 0.2114 0.0086 0.0914 -0.1890 
-0.1076 -0.6742 -0.0692 -0.6106 0.0527 -0.3847 -0.0236 


Columns 15 through 21 

-0.1854 -0.0041 -0.3831 -0.1030 -0.5743 -0.1404 -0.6633 
-0.1319 0.1225 0.2654 0.0310 0.5776 0.1159 0.8443 


Columns 22 through 28 

-0.0318 -0.7223 -0.0736 -0.7292 0.0206 -0.7679 -0.1022 
0.1240 0.9964 0.0591 1.0904 0.0730 1.1074 0.0075 


Columns 29 through 35 

-0.5967 0.0881 -0.4962 -0.0521 -0.2350 0.0872 -0.1258 
1.0717 0.0602 0.8314 -0.0597 0.5758 0.0144 0.2073 


Columns 36 through 42 

0.0442 0.0852 0.0671 0.2845 0.1206 0.5419 0.0896 
-0.1010 -0.0944 -0.0799 -0.4433 -0.1145 -0.7440 -0.1510 


Columns 43 through 49 

0.7482 0.1712 0.8393 0.0375 0.9105 0.0575 0.8713 
-1.0425 -0.1319 -1.2376 -0.1653 -1.2999 -0.0688 -1.2899 


Columns 50 through 56 

0.0530 0.8104 0.0221 0.6379 0.0741 0.3764 -0.0248 
-0.0444 -1.1347 -0.0149 -0.8797 0.0075 -0.5058 -0.0365 


Columns 57 through 61 

0.1169 -0.0111 -0.1443 -0.1061 -0.3409 
-0.0961 0.0786 0.3044 0.0666 0.7439 


Algorithm: 

The function first computes the observer/controller Markov parameters for the closed-loop 
system. Then, from the identified observer/controller Markov parameters, the individual 
system, observer, controller Markov parameters are computed and arranged in the specified 
form which is ready to be used for realization. The function is a combination of the 
functions arxc and mar_oc. 

See also: 


arxc, mar_oc, odd, ryucovar, separate 


39 


„ match 

Purpose: 

match system eigenvalues from identified forward and backward model and provide reduced 
forward and backward model in modal coordinates. 

Synopsis: 


f lamdaf,bmf,cmf,msv J,lamdab,bmb,cmb,msv_b}=match(af,bf,cf,ab,bb,cb) 

Description : 

A typical forward state space model has the form 


x(k + l) = (/4 + GC)x(k) + [B + GD 
y(k) = Cx(k) + Du(k) 



with m outputs, r inputs, and time samples dt apart. Here A, B, C, D are system matrices 
and G is referred to as the forward observer gain. The stable modes of the forward model 
are inside the unit circle, whereas the unstable modes are outside the unit circle. 

On the other hand, a typical backward state space model has the form 


x(k) = ( A ~ l + GC)x(k + 1) + [~A~'B 


GD 



«(*) 

u(k + 1) 


y(k) 


y(k) = Cx(k) + Du(k) 

The stable modes of the backward model are outside the unit circle, whereas the unstable 
modes are inside the unit circle. All matrices are the same as those described above for the 
forward model, but here G is referred to as the backward observer gain matrix. 


The matrices A, B, C, D, G, and G can be simultaneously identified by the function 
okidjb. For noise-free data, the A, B, C, D, identified either from the forward model or 
from the backward model should have an identical input-output map, implying that the 
identified system eigenvalues are identical. For noisy data, however, the system modes are 
contaminated by the noises. In addition, there are many computational (spurious) modes in 
the identified system matrices. As a result, the system matrices identified from both forward 
and backward models are somewhat different because the system modes are contaminated 
differently by noises. Note that the forward and backward models are identified by 
minimizing different residuals due to system uncertainties and measurement noises. 
Nevertheless, the system modes from both models should be reasonably close, whereas the 
computational modes may be quite different. Therefore, matching is possible to distinguish 
the system modes from the computational modes. 


Given A, B, C identified for the forward model and A~\-A~'B, C identified for the 
backward model, function match return with a reduced forward model and backward model 
in modal coordinates. The input parameters af, bf, cf represent the forward model matrices 
A, B, C, whereas ab, bb, cb represent A'\ -A~'B, Crespectively. The output vectors, 
lamdaf and lamdab, contain the matched system eigenvalues for the forward and backward 
models respectively. The matrices bmf,cmf,bmb,cmb are the corresponding input and 


40 



output matrices in modal coordinates, where the last character / means forward and b 
backward. The output vectors msv J and msvb give the msv (modal singular value, see 
function svpm) contribution to the pulse response samples for the forward and backward 
modal parameters, respectively. In addition to the eigenvalue matching, the msv 
contribution is also examined. Those modes which has higher msv contribution than the 
matched modes are also included in the reduced model. Therefore, both the reduced 
forward and backward modal models may not be the same in size. In general, the reduced 
forward modal model is larger in size for stable modes than the reduced backward modal 
model. On the other hand, the reduced backward modal model may be larger in size for 
unstable modes than the reduced forward modal model. See references for Retailed 
information. If the identified forward observer gain G and backward observer gain G are to 
be counted in the computation of the msv contribution, the input parameters bf should be 
replaced by [bfgf] and bb by [bb gb ], where gf means G and gb means G . 

Example: 

From two-input and three-output data of a three-mass-spring-dashpot system with two 
unstable modes and one stable mode, use okid and okid_b to identify a forward model and 
backward model, and then match the identified eigenvalues. The input u and output y data 
has 1 sec. sampling period, 250 data points and 10% noises. 

[af,bf,cf,df,gf]=okid(3,2,l,u,y,'batch_lq',3); 

[ab,bb,cb,db,gb]=okid_b(3,2,l,u,y,'batch_lq',3); 

ab=inv(ab);bb=-ab*bb; 

[lamdaf,bmf,cmf,msv_f,lamdab,bmb,cmb,msv_b]=niatch(af,[bf gf],cf,ab,[bb gb],cb); 

nj>=length(lamdab);n_f=length(lamdaf); 

eif=deg2hz(lamdaf,dt); 

disp([eif(:,l:2) msv_fl); 

-1.4607e-01 2.7570e-01 2.5322e-01 
-1.4607e-01 2.7570e-01 2.5322e-01 
-1.7733e-01 8.0889e-02 1.0000e+00 
-1.7733e-01 8.0889e-02 1.0000e+00 
5.5570e-01 4.4269e-01 6.9444e-02 
5.5570e-01 4.4269e-01 6.9444e-02 
lamdab=1.0 ./lamdab;bmb=-diag(lamdab)*bmb; 
eib=deg2hz(lamdab,dt); 
disp([eib(:,l:2) msv_b]); 

-1.8321e-01 2.7569e-01 2.2238e-01 
-1.8321e-01 2.7569e-01 2.2238e-01 
-1.9672e-01 8.0897e-02 1.0000e+00 
-1.9672e-01 8.0897e-02 1.0000e+00 
3.4786e-01 4.4260e-01 1.2439e-01 
3.4786e-01 4.4260e-01 1.2439e-01 

See also: 

okid, okid_b, okid_fb, svpm 

References: 

[1] Juang, J.-N. and Phan, M., "Identification of Backward Observer Markov Parameters: 
Theory and Experiments," NASA Technical Memorandum TM-107632, Langley Research 
Center, Hampton, VA., May 1992. 
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modal 

Purpose: 

Compute a reduced stable or unstable model in modal coordinates. 

Synopsis: 

[Lambda, Bm, cm\=xno<\&\{A, B ,C flag) 

Description: 

Given the discrete model 

x(k + 1) = Ax(k) + Bu(k) 
y(k) = Cx(k) + Du(k) 

with r inputs and m outputs, the equivalent model in modal coordinates is 

x»(k + \) = Ax m (k) + B m u(k) 

y = C n x(k)+Du(k) 


where A is a diagonal matrix containing the eigenvalues, A, (i‘ = l,2,...,n), of the state 
matrix A. Those eigenvalues with their length larger than 1 are known to be unstable 
modes. 

This function returns a complex vector, Lambda which contains only either stable modes 
when flag is set larger or equal to zero or unstable modes otherwise. Corresponding, the 
input matrix B is transformed to Bm and C to Cm. Note that the transmission matrix D is 
coordinate independent. This function is used in okid_b and okid_fb to distinguish 
identified unstable modes from stable modes for comparison with the identified results 
from okid. It is used in conjunction with function bk_diag to reduce an identified system 
model to a stable block diagonal form for numerical simulations to compare with real 
data. 


See also: 

okid_b, okid_fb 
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monera 

Purpose: 

Estimate variance of the ERA identified parameters. 

Synopsis: 

[eg,egm,vp,vw,mh]=mor\CT&(y,m,r,n,nm,dt,ni,nosy, 

Description: 

Monera uses the Monte Carlo approach to calculate the variance of the ERA identified 
parameters contaminated by noise. The system pulse response samples y are stored as 



’ >n(0) — 

>ml(0) 

>l r (0) 

>m,(0) ' 


>n 0) 

>ml0) 

>1,0) 

>w(l) 

> = 

>,i(2) - 

>«i(2) 

>lr(2) 

>m,(2) 


>„(/-!) - 

>m ltf-l) 

- >!,(/-!) - 

> W (/-D. 


where y 0 (t) is the /- th output at discrete time t to a unit pulse at thej'-th input. The system to 
be identified has r inputs and m outputs. The identified model order, n, is chosen by the 
user. Scalar nm specifies the row number of the Hankel matrix as described in era. Integers 
nm x m should be greater than the model order n. Scalar dt specifies the data sampling 
interval. Scalar ni specifies the number of the Monte Carlo runs used to estimate the variance 
of the ERA identified parameters. Scalar nos specifies the standard deviation of the white, 
zero-mean and Gaussian measurement noise artificially added to the pulse response samples; 
nos should be much smaller than the root mean squared value of the pulse response samples 
in y. 

[cg,e£m,vp,vw,m/i]=monera(y,m,r,rt,WM,£^m\m)s) returns the eigenvalue vector eg and 
the modal amplitude coherence vector mh of the ERA identified model from the pulse 
response samples y. The elements of vector egm are the mean values of the ni set of ERA 
identified eigenvalues from the pulse response samples with added noise (nos). Vector vp 
(vw) is the variance of the ERA identified damping (frequency) corresponding to the 
eigenvalue vector eg from ni set of Monte Carlo runs. The variable mh is the modal 
amplitude coherence. 
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Example: 

Calculate the variance of the ERA identified frequencies and dampings from the pulse 
response samples of a single input, single output second order system with sampling 
interval dt=0.2, natural frequency = 1, damping factor £ = 0.1. 

a=[0 1 1 -0.21; b=[0;l];c=[l 0];d=0; 
pt= 1 00;dt=0.2;t=[dt:dt:pt*dt 
u=zeros(pt, l);u( 1 , 1 )= 1 .0; 
y=lsim(a,b,c,d,u,t); 
rand(’normar) 

y=y+0.04*rand(pt,l); 

[eg,egm,vp,vw,r]=monera(y, 1 , 1 ,4,20,dt,50, 1 .e-5); 

eg,egm,r 

clg 

subplot(211) 

bar(vw),title('Frequency variance') 
bar(vp),title(’Damping Variance') 

Algorithm: 

The function monera uses the Monte Carlo approach [Longman89,91] and data 
correlationfJuang 88) to estimate the variance of the ERA identified frequencies and 
dampings. Depending on the number of sensors and the length of data to be used for 
system realization, monera automatically chooses either eradc or eradct to identify a system 
model with computational efficiency. 

See also: 

eradc, peradc 

References: 

[1] Longman, R. W., Bergman, M. and Juang, J. N., "Variance and Bias Confidence Criteria 
for ERA Modal Parameter Identification," Proceedings of the 1988 AAS/AIAA 
Astrodynamics Specialist Conference, Minneapolis, Minnesota, August 1988. 

[2] Longman, R. W., Lew, J. S., Tseng, D. H. and Juang, J. N., "Variance and Bias 
Computation for Improved Modal Identification Using ERA/DC," Proceedings of the 1991 
American Control Conference, Boston, MA, June 1991. 

[3] Juang, J. N., Cooper, J. E. and Wright J. R., "An Eigensystem Realization Algorithm 
Using Data Correlation (ERA/DC) for Modal Parameter Identification," Control Theory and 
Advanced Technology, Vol. 4, No. 1, pp. 5-14, 1988. 
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ocid 


Purpose: 

Identify system, observer and controller gain matrices from closed-loop test data. 
Synopsis: 

[a, b,c,d,gJ]=ocid(y,ufb,ue,p,dt, truncate, description) 

Description: 

The function ocid solves for the observer/controller Markov parameters of the following 
system 

Su/bW + u^k) 

x(k + \) = (A+GC)x(k) + {B+GD -Cl * ^ 


, 

Si 
?*- 
1 


l 

n 

L ujb(k)_ 


L-'J 


\x(k) + Du(k) 


with sampling time intervals dt apart. The data is stored as column matrices. For example, 
for a system with m outputs, the closed-loop output data matrix y contains m columns and 
as many rows as the number of data points available. The data matrix ufb contains the 
feedback control signal, ue contains the additive excitation input signal. The number p 
denotes the number of observer/controller Markov parameters to be solved for. The 
number truncate specifies the number of data points to be deleted prior to application of the 
algorithm. This value is equal to the number of time steps that is expected for the existing 
observer to converge. The description is a short descriptive tag for the current data set and 
indicates the computation procedure to be used to compute the observer and controller gain. 
If the description is set to be 'lq', it means that a, b,c,d are realized first, and then g and/ 
are computed by least-squares (lq). Any other description indicates that a, b, c, d, g, and / 
are to be realized simultaneously. The function will first return a plot of singular values for 
the user to select the desired model order. Once this is done, the function will ask whether 
or not the user wants to see plots that show the actual responses and reconstructed 
responses. Finally, a set of realized system matrices, observer gain, and controller gain 
will be returned. 

In general, if p observer/controller Markov parameters are to be solved for, then the 
maximum order of the system that can be recovered is pm , where m is the number of 
outputs. The following rules apply with regard to the number of singular values computed. 
If a flag value of zero is chosen, then the singular value plot will show (m+r)i singular 
values where i is the smallest integer such that (m+r)i is larger or equal to pm where r is 
the number of inputs. If a flag value of one is chosen, then the singular value plot will 
always show pm singular values. In either case, the user can always retain pm singular 
values to obtain a realized system that has the maximum order for a chosen value of p. 

The following information provides a better understanding of the OCID problem. Consider 
a linear discrete system of the form 



x(k + l) = Ax(k) + Bu(k) 
y(k) = Cx(k) + Du(k) 

which is operating in closed-loop. The input to the system, u(k), consists of two 
components 


u(k) = u /b (k) + u,(k) 

where Upik) denotes the feedback signal provided by an existing linear state feedback 
controller with gain F 

U/b(k) = -Fx(k) 

and u^Jk) denotes an additive excitation input for closed-loop identification. The estimated 
state x(k) is provided by an existing observer of the form 

x(k + 1) = Ax(k) + Bu(k) - G d [y(k) - y(*)] 
y(k) = Cx(k) + Du(k) 

The output y(k) is the system closed-loop response due to an excitation u e (k). OCID first 
solves for the Markov parameters 

D, an d\^^A + GC) k ' x [B + GD -G] , k = \, 2, .... p 

from which a realization of A, B, C, G, and F will be computed and returned to the user. 
Note that the matrices A, B, C, D are the system matrices and F is the existing feedback 
controller gain as described above. The matrix G , however, is another observer gain 
ass oc iated with the identified system A, B, C, D. In general, this observer matrix gain G is 
not the same as the existing observer gain G d of the closed-loop system. The matrix G is 
an observer gain for the observer given below 

x (k + 1) = Ax(k) + Bu(k) - G[y(k) - y(*)] 
y(k) - Cx(k) + Du(k) 

The function returns a, b, c, d, g, and f, which are a realization of A, B, C, D, G, and F. 

Example: 

Use test data from an aircraft flutter test. The items in italics is information prompted by the 
function ocid which has to be answered by the user. The rest is just general information 
returned by the function ocid. The following is output taken from a typical run. 

load xsamp712 

ue(l:600)=l];y(l :600)=[l;u(l :600)=[]; 

[m,n]=size(ue); 

time=dt*[0:m-l]'; 

fa,b,c,d,g,f]=ocid(y,-u,ue,30,dt,300,’ocid'); 

ERADC is used now. 

The Hankel matrix size for ERADC is 30 by 62. 


46 



SV Magnitude 


Maximum Hankel singular value = 3.213077e+02 
Minimum Hankel singular value = 2.251935e-03 



Number 


Desired Model Order (0=stop)=: 2 
Model Describes 98.9109 (%) of Test Data 
Damping(%) Freq(HZ) Mode SV MAC 
-3.2239e+00 9.1286e+00 l.OOOOe+OO 9.9827e-01 
-3.2239e+00 9.1286e+00 1.0000e+00 9.9827e-01 

Desired Model Order (0=stop)=: 4 
Model Describes 99.4891 (%) of Test Data 
Damping(%) Freq(HZ) Mode SV MAC 
2.2923e+01 1.1587e+01 5.7107e-02 9.8650e-01 
2.2923e+01 1.1587e+01 5.7107e-02 9.8650e-01 
-2.8767e+00 8.8912e+00 l.OOOOe+OO 9.9987e-01 
-2.8767e+00 8.8912e+00 l.OOOOe+OO 9.9987e-01 

Desired Model Order (0=stop)=: 8 
Model Describes 99.829 (%) of Test Data 
Damping(%) Freq(HZ) ModeSV MAC 
1.8507e+01 1.9546e+01 1.2959e-02 9.4607e-01 
1.8507e+01 1.9546e+01 1.2959e-02 9.4607e-01 
2.9880e+00 8.7176e+01 2.2222e-02 9.8971e-01 
2.9880e+00 8.7176e+01 2.2222e-02 9.8971e-01 
1.6910e+01 1.1593e+01 9.2265e-02 9.9728e-01 
1.6910e+01 1.1593e+01 9.2265e-02 9.9728e-01 
-3.3582e+00 8.8459e+00 l.OOOOe+OO 9.9995e-01 
-3.3582e+00 8.8459e+00 l.OOOOe+OO 9.9995e-01 

Desired Model Order (0=stop)=: 0 

Compare closed-loop reconst, and actual resp. (l=yes,0=no) ?:= 1 
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Act. vs. Reconst. Closed-Loop Response Pred. Open-Loop and Act Closed-Loop Resp. 


Number of sample points to reconstruct ?:= 400 
Predicted open-loop and actual closed-loop responses 



Actual vs. reconstructed closed-loop responses 
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Actual vs. reconstructed feedback control input histories 



The solid lines represent the real data whereas the dashed lines mean the reconstructed data. 
The predicted outputs are the reconstructed data from the identified system model only. 
The estimated outputs are the reconstructed data from the identified observer. It is obvious 
that the reconstructed data match the real data very well. 

Algorithm: 

Identification of the observer/controller Markov parameters for the system shown before is 
obtained using a least-squares solution. The observer/controller Markov parameters are 
identified first, from which the individual Markov parameters of the system, observer, and 
controller are computed, and used to obtain a realized state space model and the observer 
and controller gain matrices. For further details, see references. 

See also: 

arxc, mar_yoc, mar_oc, ryucovar, y_closed, separate 

Reference: 

[11 Juang, J. N. and Phan, M., "Identification of System, Observer, and Controller from 
Closed-Loop Experimental Data," Presented at the AIAA Guidance, Navigation, and 
Control Conference, Hilton Head, South Carolina, Aug. 10-12, 1992. 
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okid,okid_b,okid_fb,okid_p 

Purpose: 

Identify a state space model and its corresponding observer from test data. 

Synopsis: 


[a, ,b,c,d,g]=okid(m,r,dt,u,y, , description, p) 

[a, b,c,d,g]=okid_b(m,r,dt,y, description, p) 
{af,bf,cf,df,gf,ab,bb,cb,db,gb\=ok\dJb{m,r,dt,u,y, description#) 
[a, b,c,d,g]=okid_p(m,r,dt,y, description, p) 


Description: 

The function [A, B,C,D,G]=ok\d{m,r,dt,u,y, description, p) identifies a state space model of 
the form 

«(*)! 

m j 

y(k) = Cx(k) + Du(k) 

with m outputs, r inputs, and time samples dt apart. Here A, B, C, D are system matrices 
and G is referred to as the forward observer gain. The input variable description is a short 
descriptive tag for the current data set being analyzed and also serves as a flag which is 
described latter. The input/output time histories are stored as column matrices. For s 
experiments, the input matrix u must contain sxr columns and the output matrix y sxm. 
The number of rows in the input and output matrices equals the number of sample points. 
Initially an estimate of the number of observer Markov parameters, p, must be specified. 
For a given p, the maximum system order that can be identified is the product pxm. The 
function will prompt the user at various points for information. After computing the 
observer Markov parameters, the option to compute the prediction error is given. The 
calculation of observer Markov parameters and the output prediction error is time consuming 
when analyzing long records, therefore, it should only be used when necessary. For the 
very first run, the observer Markov parameters and related parameters including p are stored 
in the data file dokid J for function file okid, dokidb for okid_b, dokid Jb for okid_fb and 
dokid _p for okid_p. The user will be prompted if he has already stored these parameters in 
the data file. If he did, computation of these parameters will be by-passed. A plot of the 
Hankel matrix singular values is shown to aid selecting the correct system order. The 
number of non-zero singular values equals the system order. The magnitude of the Hankel 
matrix singular values, arranged in descending order, measures the state contribution. For 
noisy data, one has to make a judgement as to how many singular values to retain. After 
selecting a particular order, the percentage of the response realized by the model is computed 
using the singular values. In addition, the corresponding frequencies and damping values 
are listed with the corresponding modal amplitude coherence factors. If the model is 
acceptable, the computation is completed. The identified matrices A,B,CJ), and G are 
returned to the main program. See references for detailed information. 

The user is recommended to run the batch job for the first time so that he does not have to 
wait for the computation of observer Markov parameters which may take time for a long 
data record. The user may then come back to run the same job again interactively and use 
the existing data record for the observer Markov parameters. 


x(k + \) = (A+GC)x(k)+[B+GD -G]j 
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The user is recommended to run the batch job for the first time so that he does not have to 
wait for the computation of observer Markov parameters which may take time for a long 
data record. The user may then come back to run the same job again interactively and use 
the existing data record for the observer Markov parameters. 

okid: It simultaneously identifies A, B, C, D and G directly from input/output data. 

User's interaction is required to determine the order of the system by looking at the 
singular values plot. 

description = 'batch' ; It performs a batch job. The order of the system is 
determined internally in the era function file and thus user's interaction is not 
required. 

description = 'lq' ; It computes the system matrices. A, B, C, and D, first and then 
the observer gain G using least-squares. Researchers who are interested in 
identifying the system matrices only are recommended to use this function file. 

okid_p: a modified version of okid. It uses pulse response samples, y, to simultaneously 
identify A, B,C,D and G. No user inputs are required in this function file. 

description = 'lq' ; It computes the system matrices. A, B, C, and D, first and then 
the observer gain G using least-squares. Researchers who are interested in 
identifying the system matrices only are recommended to use this function file. 

The identified system has the following form 

x(k + l) = Ax(k) + Bu(k) 

y(k) = Cx(k) + Du(k) 

with its corresponding identified observer as 

x(k + 1) = Ax(k) + Bu(k ) - G[y(k) - y(*>] 
y(k) = C. x(k) + Du(k) 

where x (k) is the estimate of the state x(k). 

The function [A, B,C,D l G]=ckid_b(m,r,dt,u,y, description^) identifies a state space model 
of the form 


x(k) = (A~' +GC)x(k + \) + [-A~'B GD 


u(k) 



u(k + 1) 


y(k) 


y(k) = Cx(k) + Du(k) 


with m outputs, r inputs, and time samples dt apart. All the input and output parameters of 
this function are the same as jhose described above for the forward observer. Here A, B, C, 
D are system matrices and G is referred to as the backward observer gain matrix. See 
references for detailed information. 


okid_b: It simultaneously identifies A, B, C, D and G directly from input/output data. 
User's interaction is required to determine the order of the system by looking at the 
singular values plot. 

description = 'batch' ; It performs a batch job. The order of the system is 
determined internally in the era function file and thus user's interaction is not 
required. 

description = 'lq' ; It computes the system matrices, A, B, C, and D, first and 
then the observer gain G using least- squares. Researchers who are interested in 
identifying the system matrices only are recommended to use this function file. 

The identified system is the same as above, whereas the observer becomes 

x(k) = A~ l x(k + 1) + A~ l Bu(k ) - G[y(k + 1) - y(k + 1)] 

y(k) = Cx(k) + Du(k) 

where x(k) is the estimate of the state x(k). Note that the backward approach identifies the 
inverse of the system state matrix whereas the forward approach identifies the system state 
matrix directly. The advantage of the backward approach is that all the identified spurious 
modes tend to be stable and the strong system modes to be unstable. Therefore when the 
identified state matrix is inverted, the strong system modes become stable but the 
computational modes become unstable. Nevertheless, experiences suggest that the 
identified results are somewhat underestimated particularly when noises are high. To do 
numerical simulations, functions modal and bk_diag may be used to reduce the identified 
model to a stable model in the real domain. 

The function [Af ,Bf,Cf,Df,Gf,Ab ,Bb ,Cb,Db ,Gb]=okid_fb(m,r,dt,u,y , description#) 
identifies two state space models simultaneously using both the forward and backward 
approach. The small cases / and b behind the capital characters means forward and 
backward respectively. In this function, comparison of stable system modes in terms of 
frequencies and dampings from these two models is provided to the user for his judgement 
on how accurate the system modes are. The comparison is based on the error between the 
forward system eigenvalues and backward system eigenvalues, and their singular value 
contributions to the pulse response samples. The comparison provides the reduced output 
forward model [Af,Bf,Cf,Df,Gf\ and backward model [Ab,Bb,Cb,Db,Gb] where the state 
matrices ,4/ and Ab are both in block-diagonal form (see function bk_diag). Note that the 
reduced modal observer may not be stable, since the modal reduction is not optimal in 
general. In this case, the user is recommended to use function okid to identify a system 
model and then reduced by other methods such as the optimal projection or balanced 
coordinates. The function okid_fb is strongly recommeded for those users who do not care 
about the observer identification. With the same order of observer markov parameters, it is 
believed that the forward approach provides the identified results with better accuracy. The 
forward model is in general larger in size than the backward model. The backward 
approach may miss some system modes particularly with light damping. However, it 
provides information of strong system modes. 
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i. Smoh. Error Norm. Pred. Error 


Example : 


Use test data from a truss structure. The items in italics is information prompted by the 
function which has to be answered by the user. The rest is just general information returned 
by the function. The following is output taken from a typical run. 


load x sample 

la,b,c,d,m]=okid_fb(n,r,dt,u,y, , okid_fb',20); 

Total number of sample points = 2000 
Number of experiments in file = 1 
Number of inputs = 2 
Number of outputs = 2 

Compute Observer Paramters For Data Set Number 1 
Time (min) to compute parameters 4.1 14 

Have you run OKID_FB with the same data & P before (l=yes,0=no) ?:= 0 

Compute Forward and Backward Error (l=yes,0=no)? =: 1 
Compute Prediction Error For Data Set Number 1 
Forward Square Fitting Error Normalized 
7.263 le-02 -4.8648e-02 
-4.8648e-02 1.121 le-01 
Backward Square Fitting Error Normalized 
1.5161e-01 -4.6554e-02 
-4.6554e-02 1.4184e-01 
Time (min) to Markov parameters 0.1517 


0.02 



. 0 . 02 1 * * * — * * ■ ■ ■ * 1 

0 200 400 600 800 1000 1200 1400 1600 1800 2000 

Time Steps 


£ 


0.02 



-0 02 ‘ * i i i i i : . . . l 

0 200 400 600 800 1000 1200 1400 1600 1800 2000 


Time Steps 


53 





SV Magnitude 


THE FOLLOWING COMPUTES A DISCRETE MODEL FROM A FORWARD OBSERVER. 



Number 


Desired Model Order (0=stop)=: 26 
Model Describes 99.9924 (%) of Test Daia 
Damping(%) Freq(HZ) Mode SV MAC 
9.8448c+00 9.7606C+01 3.0499c-03 9.9677e-01 
9.8448e+00 9.7606e+01 3.0499c-03 9.9677e-01 
7.4294c+00 l.l862e+02 2.2901c-02 9.9957e-01 
7.4294e+00 1.1862e+02 2.2901e-02 9.9957e-01 
5.1193c+00 1.0573e+02 7.2535c-03 9.9891e-01 
5.1193c+00 1.0573c+02 7.2535e-03 9.9891c-01 
1.7537e+01 2.2577c+01 2.3866c-03 9.9869e-01 
1.7537e+01 2.2577c+01 2.3866c-03 9.9869e-01 
1.1365e+01 3.3578c+01 9.4557e-04 9.8801e-01 
1.1365e+01 3.3578c+01 9.4557c-04 9.8801e-01 
1.7740c+01 1.9096c+01 1.4288c -03 9.9687e-01 
1.7740c+01 1.9096e+01 1.4288e-03 9.9687e-01 
1 .9292e+00 1.1410e+02 4.451 le-02 9.9999e-01 
1.9292e+00 1.1410e+02 4.451 lc-02 9.9999e-01 
3.4358e+00 6.1287e+01 1.6785c-03 9.9893e-01 
3.4358C+00 6.1287e+01 1.6785e-03 9.9893c-01 
1.4563e+00 4.659 lc+01 1.2229c-02 9.9996e-01 
1.4563e+00 4.6591c+01 1.2229c-02 9.9996e-01 
6.7864e-01 7.424 lc+01 4.4470c-03 9.9929e-01 
6.7864c-01 7.424 lc+01 4.4470c-03 9.9929e-01 
4.2429c-01 4.8646e+01 6.9666c-02 1.0000e+00 
4.2429c-01 4.8646c+01 6.9666e-02 1.0000e+00 
3.1829c-01 7.2733e+00 4.8384c-01 1.0000e+00 
3.1829e-01 7.2733c+00 4.8384e-01 1.0000c+00 
3.9005c-01 5.8479e+00 1.0000c+00 1.0000c+00 
3.9005C-01 5.8479c+00 l.()000c+00 1.0000c+00 

Desired Model Order (0=stop)=: 0 
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SV Magnitude 


THE FOLLOWING COMPUTES A DISCRETE MODEL FROM A BACKWARD OBSERVER. 


Hankel Matrix Singular Values 



Number 


Desired Model Order (0=stop)=: 26 
Model Describes 100 (%) of Test Data 
Damping(%) Frcq(HZ) Mode SV MAC 
8.7805e+00 1.0444e+02 1.5233e-02 9.9862e-01 
8.7805C+00 1.0444e+02 1.5233e-02 9.9862e-01 
6.9248e+00 1.1904e+02 5.8056e-02 9.9986e-01 
6.9248e+00 1.1904e+02 5.8056e-02 9.9986e-01 
5.7392e+00 1.0529e+02 2.2203e-02 9.9896e-01 
5.7392e+00 1.0529e+02 2.2203c-02 9.9896e-01 
1.9648c+00 1.1406c+02 8.8767c-02 1.0000e+00 
1.9648C+00 1.1406c+02 8.8767c-02 1.0000e+00 
7.0426e+00 3.1259e+01 2.6690e-03 9.8609e-01 
7.0426e+00 3.1259e+01 2.6690e-03 9.8609c-01 
2.6303c+00 6.1469e+01 2.4191e-03 9.9780e-01 
2.6303e+00 6.1469c+01 2.4191e-03 9.9780e-01 
4.7085c+00 1.9505c+01 1.6134e-02 9.9995e-01 
4.7085e+00 1.9505c+01 1.6134e-02 9.9995c-01 
1.8700C+00 4.6542e+01 2.1705c-02 9.9976e-01 
1.8700c+00 4.6542e+01 2.1705c-02 9.9976c-01 
6.9673e-01 7.4489c+01 4.6185c-03 9.9365e-01 
6.9673c-01 7.4489e+01 4.6185e-03 9.9365e-01 
4.2593c-01 4.8600e+01 8.9072c-02 9.9999e-01 
4.2593c-01 4.8600c+01 8.9072c-02 9.9999e-01 
-8.8086e-01 5.8583c+00 1.0000c+00 1.0000e+00 
-8.8086c-01 5.8583e+00 1.0000e+00 1.0000e+00 
-9.381 lc-01 7.3066c+00 3.9005e-01 9.9999e-01 
-9.381 le-01 7.3066e+00 3.9005e-01 9.9999e-01 
-1.9640e+01 1.4194e+01 3.1534e-02 1.0000e+00 
-1.9640e+01 1.4194c+01 3.1534e-02 1.0000e+00 

Desired Model Order (0=stop)=: 0 
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COMPARISON OF FORWARD AND BACKWARD IDENTIFICATION 


Forward Identification 


Damping(%) 

1.7740e+01 

1.7740e+01 

3.4358e+00 

3.4358e+00 

1.7537e+01 

1.7537e+01 

1.4563e+00 

1.4563e+00 

6.7864e-01 

6.7864e-01 

9.8448e+00 

9.8448e+00 

5.1 193e+00 

5.1193e+00 

4.2429e-01 

4.2429e-01 

7.4294e+00 

7.4294e+00 

1.9292e+00 

1.9292e+00 

3.1829e-01 

3.1829e-01 

3.9005e-01 

3.9005e-01 


Freq(hz) 
1.9096e+01 
1.9096e+01 
6.1287e+01 
6.1287e+01 
2.2577e+01 
2.2577e+01 
4.659 le+01 
4.659 le+01 
7.424 le+01 
7.424 le+01 
9.7606e+01 
9.7606e+01 
1.0573e+02 
1.0573e+02 
4.8646e+01 
4.8646e+01 
1 . 1 862e+02 
1 . 1 862e+02 
1.141 0e+02 
1.1410e+02 
7.2733e+00 
7.2733e+00 
5.8479e+00 
5.8479e+00 


Mode SV 
1.8532e-03 
1.8532e-03 
3.2029e-03 
3.2029e-03 
3.2734e-03 
3.2734e-03 
1.3625e-02 
1.3625e-02 
1.7323e-02 
1.7323e-02 
2.8269e-02 
2.8269e-02 
5.7957e-02 
5.7957e-02 
7.3772e-02 
7.3772e-02 
2.0175e-01 
2.0175e-01 
3.2723e-01 
3.2723e-01 
4.4107e-01 
4.4107e-01 
1.0000e+00 
1.0000e+00 


Backward Identification 
Damping(%) Freq(hz) Mode SV 


1.9640e+01 
1.9640e+01 
9.381 le-01 
9.381 le-01 
8.8086e-01 
8.8086e-01 


1.4194e+01 

1.4194e+01 

7.3066e+00 

7.3066e+00 

5.8583e+00 

5.8583e+00 


3.208 le-02 

3.208 le-02 

3.4642e-01 

3.4642e-01 

1.0000e+00 

1.0000e+00 


These are the modes selected from for wa r d and backward models 
You should also examine the list of modes from the forward model 
to see if there are some other modes left out from the above table. 
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DATA RECONSTRUCTION FROM THE IDENTIFIED FORWARD MODEL 
Compare Recons. Output and True output (l=yes,0= no) ?:= 1 


Number of Sample Points to Reconstruct ?:= 1000 

Comparison For Data Set Number No. 1 

The following figures show predicted and real outputs 
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Predicted output Predicted output 


DATA RECONSTRUCTION FROM THE IDENTIFIED BACKWARD MODEL 
Compare Recons. Output and True output (l=yes,0=no) ?:- 1 

Number of Sample Points to Reconstruct ?:= 1000 

Comparison For Data Set No. 1 

The following figures show smoothed and real outputs 
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The solid lines in the above figures represent the real data whereas the dashed lines mean the 
reconstructed data. The predicted outputs are the reconstructed data from the identified 
system model only. The estimated outputs are the reconstructed data from the identified 
observer. It is obvious that the reconstructed data from the identified forward observer 
match the real data much better than that from the identified model only. With the same 
order chosen which is 26 in this example, the forward results are somewhat better than the 
backward results (see the reconstructed predicted outputs). In this example, only three 
stable modes (after the identified state matrix was inverted) are identified from the backward 
approach. 


Algorithm: 

Identification of the pulse response (Markov parameters) for the observer system shown 
before is obtained using singular value decomposition. The forward or backward observer 
is identified first, as opposed to the system itself, such that the observer has all its poles 
placed at the origin. From this, the system Markov parameters (pulse response) are 
recovered and used in system realization. Theoretically, there are only a specific number of 
independent system Markov parameters for a finite set of observer Markov parameters. 
Therefore a minimum number of system Markov parameters may be used in the function era 
or eradc to minimize the computational time in identifying the system. Nevertheless, it 
seems from experience that a little larger number than the minimum one for the system 
Markov parameters help a little bit for the identification of system with very low damping. 
The realization algorithm provides a state space model and permits the evaluation of different 
system orders. Order selection is guided by the singular values of a Hankel matrix, but for 
test data it is up to the user to decide. For details see references. 

See also: 

arx_b, arx_bat, arx_fb, era, eradc, k_abcd, mar_com, match, pred_err, svmp, uy_stack, 
yucovar 
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p2m 

Purpose: 

Rearrange pulse response time histories in the form of Markov parameters sequence. 

Synopsis: 

y o =p2m(V) 

Description: 

Given a discrete model 

x(k + \) = Ax(k) + Bu(k) 

y(k) = Cx(E) + Du(k) 

with r inputs and m outputs, the pulse response samples are typically stored as 


y p =[y x y* ••• y~] 

where y, = response samples due to a unit pulse at the i-th input. Each y i has m columns 
and t rows where l is the length of data. The sequence of system Markov parameters is 
defined as 


Y 0 =[D CB CAB ••• ca 12 b] 

The sequences Y 0 and Y p are equivalent in the sense that both represent pulse response 
samples. The function converts Y p to Y 0 . Note that Y p is the sequence used in the eradc 
and era functions. 

See also: 

m2p, okid, okid_b, okid_fb 
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peradc 


Purpose: 

Calculate variance and bias of the ERA/DC identified parameters. 

Synopsis: 

[eg,vp,vw,bp,bw,sg,vsg,bsg]=peTadc(y,n,n m,nr,dt)\ 

Description: 

The function peradc calculates the variance and bias of the single-input and single-output 
eradc identified parameters from pulse response samples. The column vector y is the pulse 
response samples in which the t'th element of y is the unit pulse response at discrete time i. 
Scalar n specifies the identified model order. Scalar nm(nr ) specifies the number of rows 
(columns) of a Hankel matrix. In peradc, nr is required to be larger than or equal to nm. 
Scalar dt specifies the data sampling interval. 

[. eg ,vp ,vw ,bp ,bw ,sg,vsg ,bsg]=pzradc(y ,n,nm,nr ,dt ) returns the variance and quadratic 
bias of the eradc identified parameters. The elements of vector eg are the eigenvalues of the 
eradc identified model. Vector vp (vw) is the variance of the eradc identified dampings 
(frequencies) corresponding to eigenvalue vector eg. Vector bp ( bw ) is the quadratic bias of 
the eradc identified dampings (frequencies) corresponding to eigenvalue vector eg. Vector 
sg is the singular value vector of the correlation matrix in eradc. The elements in vector vsg 
( bsg ) is the normalized variance (quadratic bias) of the singular values in sg. 

Example: 

Calculate the variance and bias of the eradc identified parameters from the noisy pulse 
response samples of a single-input single-output second-order system with sampling 
interval dt=0.2, natural frequency ft), = 1, and damping factor £ = 0.1. 


a=[0 1;-1 -0.2]; b=[0;l];c=[l 0];d=0; 
pt=100;dt=0.2;t=[dt:dt:pt*dt]'; 
u=zeros(pt, 1 );u( 1 , 1 )= 1 .0; 
y=lsim(a,b,c,d,u,t); 
rand('normal') 
y=y+0.04*rand(pt, 1 ); 

feg,vp,vw,bp,bw,sg,vsg,bsg]=peradc(y,4,10,20,dt); 

eg 

clg 

subplot(221) 

bar(vw),titleCFrequency variance') 
bar(vp),title('Damping variance’) 
bar(bw),titleCFrequency bias') 
bar(bp),title(’Damping bias') 
pause 
clg 

subplot(121) 

bar(vsg),title(' Variance of singular values') 
bar(bsg),title('Quadratic bias of singular values') 
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Algorithm: 

The function peradc uses the algorithm from |Longman89,91]. In peradc the perturbation 
theory is applied to estimate the variance and quadratic bias of the ERA/DC identified 
parameters. 

See also: 

monera, e rad c 

References: 

[1] Longman, R. W., Bergman, M. and Juang, J. N., "Variance and Bias Confidence Criteria 
for ERA Modal Parameter Identification," Proceedings of the 1988 AAS/AIAA 
Astrodynamics Specialist Conference , Minneapolis, Minnesota, August 1988. 


[21 Longman, R. W., Lew, J. S., Tseng, D. H. and Juang, J. N., "Variance and Bias 
Computation for Improved Modal Identification Using ERA/DC," Proceedings of the 1991 
American Control Conference, Boston, MA, June 1991. 
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pred efb, prederb, pred_err 

Purpose: 

Compute the prediction error based on identified parameters. 

Synopsis: 

[errorf,yhatf]=pred_en(u,y, Yf,pjlag) 

[errorb,yhatb]=pTcd_err(u,y,Yb,pflag) 

[errorf.yhatf, errorb,yhatb]=pred_cn(u,y,Yf,Yb,pJlag) 

Description: 

Given the input matrix u of dimension l x r, the output matrix y dimension / x m, and an 
estimate of p forward observer Markov parameters stored in Yf, the forward prediction error 
computed in matrix form is 


where 



y 7 = ly(0)y(i)-y(/-i)] 

Y,=\D CB CAB — CA pl B] 

'k(O) m(1) u( 2) -. u{l- 1) 

v(0) v(l) v(l- 2) 

V ~ : : 

v(0) v(l-p- 1) 


K0 = 


u(i) 

yd) 




,/-i 


See functions okid and arx_bat for the definition of matrices shown above. The number of 
samples is / and the system has m outputs and r inputs. The computation is performed 
within a loop to reduce storage requirements at the expense of computation time. For cases 
where memory is not a problem flag=\ yields faster computation time. The function 
pred_err returns the prediction error in the vector errorf and the estimated measurement y in 
the vector yhatf. 

Given an estimate of p backward observer Markov parameters stored in Y b , the backward 
smoothing error computed in matrix form is 

e> = y b ~Y b V b 

where 
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y b =[y(0)y(\>-y(l-p-2)) 

Y„=[D + CB CA^GD -G] 

r«(0) k(1) u( 2 ) 

„ v(p) v(p + 1) v(p + 2) 


CA'~ 2 [AB + GD -G] 

••• u(l- p- 2) 

— v(/-2) 


KO = 


v(l) 

■«(/)! 

_y(Oj 


v(2) 


1 = 0 , 1 , 


v(3) 


,/-l 


v(/-p-l) 


••• c[ab + gd -G]] 



See functions arx_b, arx_fb, okid_b and okid_fb for definition of the matrices shown 
above. The function pred_erb returns the backward smoothing error in the vector errorb 
and the smoothed measurement y in the vector yhatf. 


Due to the strong similarity between the matrices V f and V b , the forward and backward 
errors may be computed simultaneously. The function pned_efb returns the prediction error 
in the vector error/, the smoothing error in errorb, the estimated measurement y in yhatf, 
and the smoothed measurement y in yhatb. 

Example: 

m= 1 ;r= 1 ;p=2;L=5;flag= 1 ; 

a=[0 -0.16; 1 -1]; b=[0 11’; c=[0 1 J; d=0; G=[0.16 1]’; 

u=rand(L,r); 

y=dlsim(a,b,c,d,u); 

abar=a+G*c;bbar=fb+G*d -G]; 

Y=[d c*bbarc*abar*bbar]; 

[error,yhat]=pred_err(u,y,Y,p,flag); 

See also: 

arx_bat 
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Purpose: 


pulse 


Compute pulse response histories from general input and output data. 

Synopsis: 

[ys,yo]=puhe(m,r,dt,u,y,p,n _pulse, description); 

Description: 

The function pulse computes the unit pulse response samples (Markov parameters) from 
input and output time histories by using a time domain approach. The system input and 
output histories u(t) and y(t) must be stored as follows 



«n(0) 

u r i(0) 

Kli(O) 

U rs ( 0) ' 


MO 

MO 

MO 

Ml) 

u = 

m u (2) 

U r 1<2) 

m 1j(2) 

M2) 


.«u(/-1) - 

un(l- 1) 

"■ u u (l-V) 

Ursd- 1). 


yn(0) 

>W(0) 

yi*(0) 

yms( o) 


>-ii 0) 

y* id) 

yu <0 

^(1) 

y = 

y„(2) — 

y m i(2) 

^u(2) 

yms(2) 


yn(/-l) - 

JW0-D 

yii(f-i) ■■■ 

ymsV- 1) 


where r is the number of inputs and m is the number of outputs, and u^(t)(yJt)) is the i-th 
input (output) of the y'-th test at discrete time t. Multiple experiments are allowed in using 
this function file. The input variable p is the desired number of independent observer 
Markov parameters to be identified from input and output time histories. Given the desired 
number p, the maximum number of the system order is p*m. The input variable n _pulse is 
the desired length of the pulse response time histories to be computed. The input variable 
description is a short descriptive tag for the current data set being analyzed and also serves 
as a flag. If description = 'inverse’, it computes the pulse response histories for the 
backward model which can be used to realize the inverse of a system state matrix (see 
description for function okid_b). This function works for stable and unstable systems. All 
calculations are performed in the time-domain. The output of this function is the system 
pulse response histories ys of dimension n jjulse by mr and the observer pulse response 
histories yo. of dimension n jjulse by mm. 

Example: 

This example is to identify the pulse response time histories of a three-mass-spring-dashpot 
system from two-input and one-output data. 


kl=1.0;k2=2.0;k3=3.0; 

ml=l .0;m2=1.0;m3=l .0;ratio=2*0.005; 

K=[kl+k2 -k2 0 ; -k2 k2+k3 -k3; 0 -k3 k3]; 
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Khalf=sqrtm(K);Damp=ratio*Khalf; 

Ac=[zeros(3,3) eye(3,3);-K -Damp]; 

Bc=[zeros(3,2);l 0; 0 1; 0 0]; 

C=[ zeros (1,5) 1]; 
dt=l .0;pt=200;p=10; 

[m,n]=size(C);[n,r]=size(Bc); 

D=zeros(m,r); 

t=[dt:dt:pt*dt]'; 

[A,B]=c2d(Ac,Bc,dt); 

rand('normal'); 

u=rand(pt,r);y=dlsim(A,B,C,D,u); 
[ys,yo]=pulse(m,r,dt,u,y,p,pt, 'pulse'); 

%[ys,yo]=pulse(m,r,dt,u,y,p,pt,’inverse');%for backward only 


Algorithm: 

This is a time-domain approach. First, the forward or backward observer Markov 
parameters are identified (see arx_bat, arx_b), as opposed to the system itself, such that the 
corresponding observer has all its poles placed at the origin. From these identified observer 
Markov parameters, the system pulse response is recovered (see mar_com). For details see 
references. 

See also: 

arx_bat,arx_b,mar_com,m2p,okid,okid_b 

References: 

]1] Chen, C. W., Huang, J.-K., Phan, M„ and Juang, J.-N.,"Integrated System Identification 
and Modal State Estimation for Control of Large Flexible Space Structures," Journal of 
Guidance, Control and Dynamics, Vol. 15, No. 1, Jan.-Feb. 1992, pp. 88-95. 

[2] Phan, M., Horta, L. G., Juang, J.-N., and Longman, R. W„ "Linear System Identification 
Via an Asymptotically Stable Observer," Proceedings of the AIAA Guidance, Navigation and 
Control Conference , New Orleans, Louisiana, Aug. 1991, pp. 1180-1194, and NASA 
Technical Paper 3164, 1991, and to appear in the Journal of Optimization Theory and 

Application. 


[3] Juang, J.-N., Phan, M., Horta, L. G.,and Longman, L. G., "Identification of Observer and 
Kalman Filter Markov Parameters: Theory and Experiments," Proceedings of the AIAA 
Guidance, Navigation and Control Conference, New Orleans, Louisiana, Aug. 1991, pp. 
1195-1207. 

[4] Phan, M., Juang, J.-N., and Longman, R. W., "On Markov Parameters in System 
Identification," NASA Technical Memorandum TM-104156, Langley Research Center, 
Hampton, VA., Oct. 1991. 

[5] Juang, J.-N. and Phan, M., "Identification of Backward Observer Markov Parameters: 
Theory and Experiments," NASA Technical Memorandum TM-107632, Langley Research 
Center, Hampton, VA., May 1992. 
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ryucovar 


Purpose: 

Compute the left correlation matrix associated with the feedback control input for 
observer/controller identification. 

Synopsis: 

[ufbVt] = ryuco\ar(y,ufb,ue,p) 

Description: 

The data is stored as column matrices. For a system with m outputs, the closed-loop 
output data matrix y contains m columns and as many rows as the number of data points 
available. The data matrix ufb contains the feedback control signal, ue contains the additive 
excitation input signal in the same format as y. Let the given data be arranged as 

T = [y(0) y( 1) y(2) - y(N- 1)] 

^=[V°) V 1 ) M 2 > - V*- 1 )] 

/? = [«.((>) u t ( 1) u t ( 2) - u e (N- 1)] 

The function ryucovar computes the following correlation matrix 

ufbVt = U fl,V T 


where the data matrix V is defined as 



m(0) m(1) u( 2) 

U(p ) u(p + 1) 

- u(N - 1) ' 


z(0) z(l) 

— z{p- 1) z(p) 


- z(N- 2) 

V = 

z(0) 

z(p- 2) zip- 

• ♦ • 

• V • 

-1) 

— z(N - 3) 




• • 
2(0) Z(l) 

- z(N-p-l) 

and z(k) is defined as 







z(k) = 

’UflW + U'ik) 

y«) 




As indicated, the number N specifies the number of data points to be used in the 
computation. N can be less than the total number of data points available in ufb, ue, and y. 
The matrix product returned by the function is used in the computation of the 
observer/controller Markov parameters in the function odd. 
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Example: 


ufb=[l -2 1 3 1 345 6]'; 
y=[0 3 2 1 3 1 4 4 61'; 
ue=[-l 0 2-1 2 3 507]’; 
[ufbVt] = ryucovar(y,ufb,ue,2) 
ufbVt = 

70 42 23 15 26 

Algorithm: 


To save memory space, the summations involved in the product (/^V r are performed using 
inner matrix product multiplication. 

See also: 

ocid, arxc, mar_yoc, mar_oc, y_closed, separate 
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separate 

Purpose: 

Separate a given matrix sequence into two sequences with prescribed formats. 

Synopsis: 

[ Y 1 ,Y2] =separate(F,<?,m 1 ,m2) 

Description: 

Given a matrix Y in the following format 

Y = [Y { ( 1) F,(2) Y 2 ( 1) Y 2 ( 2) Y„( 1) Y n ( 2)] 

this function returns FI and F2 which are 

n=[K,( 1) r 2 (i) ... n(D] 

K2 = [K,(2) F 2 (2) -. F*(2)] 

Each Ji(l) has dimensions <7 x ml , and each Y t (2) has dimensions q x m2 . 

Example: 

Y=rand(2,9) 

Y = 

Columns 1 through 7 

0.9304 0.5269 0.6539 0.7012 0.7622 0.0475 0.3282 

0.8462 0.0920 0.4160 0.9103 0.2625 0.7361 0.6326 

Columns 8 through 9 

0.7564 0.3653 

0.9910 0.2470 

[Y1 ,Y2]=separate(Y,2,l ,2) 

Y1 = 

0.9304 0.7012 0.3282 

0.8462 0.9103 0.6326 

Y2 = 

0.5269 0.6539 0.7622 0.0475 0.7564 0.3653 

0.0920 0.4160 0.2625 0.7361 0.9910 0.2470 

See also: 

arxc, mar_yoc, mar_oc, ocid, ryucovar, y_closed 
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svpm 

Purpose: 

Compute modal observability matrix and singular values of the modal participation to the 
pulse response samples. 

Synopsis: 

[svm, obsm]=svpm(lambda,bm,cm,n) 

Description: 

Consider the discrete model in the modal coordinates 

x m (k + 1) = A x m (k) + B m u(k) 

y = C m x(k) + Du(k) 


with r inputs and m outputs, where A is a diagonal matrix containing the eigenvalues, 
%.(i = 1,2, of the system matrix. The modal observability matrix is computed by 


obsm = 




Function [svm, obsm]=s\pm(lambda,bm,cm,n) returns the complex modal 
observability matrix and a normalized singular value vector svm (see algorithm) to 
quantify the importance of each individual mode, for given system eigenvalues in the 
vector lambda, the modal input matrix bm, modal output matrix cm and the desired 
length n. Note that the maximum singular value is used to normalize the vector svm. 

Example: 

A three-mass-spring-dashpot system from two-input and three-output data is used. 


kl = 1 .0;k2=2.0;k3=3.0; 
m 1 = 1 .0;m2= 1 .0;m3=l .0;ratio=2*0.005; 

K=[kl+k2 -k2 0 ; -k2 k2+k3 -k3; 0 -k3 k3]; 
Khalf=sqrtm(K);Damp=ratio*Khalf; 

Ac=[zeros(3,3) eye(3,3);-K -Damp]; 

Bc=fzeros(3,2);l 0; 0 1; 00]; 

C=[zeros(3,3) eye(3,3)]; 

dt=1.0;pt=60; 

lA,B]=c2d(Ac,Bc,dt); 

[V,lambda]=eig(A); 

bm= V\B;cm=C*V ;lambda=diag(lambda); 

[sv,obsvm] = svpm(lambda,bm,cm,6); 
sv' 

sv’ = 0.2898 0.2898 1.0000 1.0000 0.1506 0.1506 
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Algorithm: 


For a linear system, the map from input u to output y can be fully described by the 
Markov parameter sequence 

Y.=[D C,B. C.AB. - C.A'-’B.] 

This sequence is coordinate independent and unique. Let the input and output matrices 
be partitioned as 



where n is the number of modal coordinates, /?/ (i-1 a row vector of length r 
and c/ a column vector of length m. Each individual Markov parameter can then be 
written as a combination of n components contributed from different modal coordinates. 

C.AB. = 

1=1 

Therefore, each coordinate has a sequence of Markov parameters described as follows 

Y* =[0 CM CM - CX 2 B,}, i = 1,2 /I 

The total Markov parameter sequence becomes 

y. 0 =[£> 0 0 0] 

i=0 

From this representation, it is obvious that each modal coordinate contributes to the 
pulse response sample by the individual modal sequence Y^, which can be quantified 
by taking its maximum singular value, i.e., 

svm = ^|^\(l+\X i \+\X 2 i \+...+\2. t ~ 2 \)^J\b~\ * 

where \Xj\ is assumed to be less than 1. 

See also: 

era, eradc, match 
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svra 

Purpose: 

Identify a state space model from input and output datavia a state vector realization 
algorithm. 

Synopsis: 

| a,ft,c,</,eg,5g]=svra(tt,y,«,rtw,/ir,rfr,A:); . . , r 

Description: 

The function svra identifies a state-space model of a multi-input and multi-output linear, 
time- invariant system from a set of rich input response data. The i-th row vector of matrix u 
( y ) is the system input (output) at discrete time i. Scalar n specifies the identified model 
order. Scalar nm specifies the number of sample shift in constructing the rows of a Hankel 
matrix and it is required to be even, whereas nr specifies the column number of the Hankel 
matrix. Integers nmxm(m is the number of outputs) and nr are required to be not smaller 
than the chosen model order n. Also, n + nm x r(r is the number of inputs) needs to be 
smaller than nr. Scalar dt denotes the data sampling interval, and the data used for realization 
starts from the k- th discrete time. [A ,B ,C ,D ,eg ,sg]=svra(u,y,n,nm,nr,dt,kflag) returns an 
nth-order linear, time-invariant identified discrete model: 

*(* + !)= Ax(k) + Bu(k) 

y(k) = Cx(k) + Du(k) 

Vector eg contains modal parameters of the identified model including frequencies (Hz) in 
the first column and damping ratios (%) in the second column. The third column of eg gives 
the eigenvalues of the corresponding continuous time model. Note the identified discrete 
model can be easily transformed to a continuous time model. Vector sg, whose elements are 
singular values of the Hankel matrix, can be used as reference to choose the model order n. 

Examples: 

i) Use svra to identify a system model from a set of random input response data for a single 
input, single output second order system with sampling interval dt= 0.3, natural frequency 
co H = 1 , and damping factor £ = 0.1. 

ii) Transfer the identified discrete model to a continuous-time model and plot the original 
system output and the output from the svra identified model. 



a=[0 1;-1 -0.2]; b=[0;l];c=[l 0];d=0; 

pt=100;dt=0.4;t=[dt:dt:dt*pt]'; 

rand('normal'); 

u=rand(pt,l); 

y=lsim(a,b,c,d,u,t); 

[al,bl,cl ,dl ,sg,eg]=svra(u,y,2,20,40,dt,10); 
y 1 =dlsim(a 1 ,b 1 ,c 1 ,d 1 ,u); 
y2=[y yll; 
clg 

plot(y2),title(’Randoni input response') 
pause 

error=y-yl; 

clg 

plot(error),title('Error’) 

pause 

sg 


Algorithm: 

The function svra uses the state vector realization algorithm from [Moonen89 Lew91]. It 
uses output and input data to form the measurement matrix as 


Uk 

u k + 1 

u k+nr-l 

yk 

1 

yk+nr-\ 

1 

u k + 2 

U k+ nr 

y*+ 1 

y*+ 2 

yk+nr 

* • 
m # 

U k+run - 1 

^k+nm 

Mk+nm+nr-2 

^k+nm-l 

yk+tun 

*" yk+nm+nr-2_ 


where Ui and yj denote r-dimensional input vector and m-dimensional output vector at time 
respectively. In svra, the SVD of H is used for a state vector realization [Moonen89] , i.e. 
Xk+i. x k+M.-' x *+M. x k+j' J >>L It then uses the following equation 

x **hA 

••• J l c D k u ^ “*+>-> J 

to identify a state-space model. 

References: 

[1] Moonen, M., Demoor, B., Vandenberghe, L. and Vandewalle, J., "On- and Off-Line 
Identification of Linear State-Space Models," International Journal of Control, Vol. 49, 
1989, pp 219-232. 

[2] Lew, J. S., Juang, J. N. and Longman, R. W., "Comparison of Several System 
Identification Methods for Flexible Structures," Proceedings of the 32nd Structures, 
Structural Dynamics and Materials Conference, Baltimore, MD, April 1991, pp. 2304- 
2318. 
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uystack 

Purpose: 

Compute a stacked matrix with inputs and outputs. 

Synopsis: 

[v.v/]=uy_stack(«,y,p) 

Description: 

Given a set of data with l samples, r inputs, m outputs, and an assumed system order of 
p*m, the function stacks a matrix using the inputs and outputs as follows 


V = 


W 

v(/)= t = 0,1,..., 

bO)J 


U( 2) - 

U{1- 1) " 

v(l) ••• 

v(/-2) 

v(0) ■” 

v(l-p- 1) 

0 

1 



The matrix dimension is [O' + m )P + r]xl 


Example: 


u=[0 1 2 3 4 5]’; 
y=[6 7 8 9 10 11]'; 

P=2; 

[vst]=uy_stack(u,y,p) 
vst = 

0 1 2 3 4 5 

0 0 1 2 3 4 

0 6 7 8 9 10 

0 0 0 1 2 3 

0 0 6 7 8 9 


See also: 


arx_bat 
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y_cIosed 

Purpose: 

Reconstruct feedback control input and closed-loop response using odd-identified system, 
observer gain, and controller gain matrices. 

Synopsis: 

[ufb_rec,y_rec]=y_closed(A,B,C,D,G,F,y,ufb,ue) 

Description: 

The function reconstructs the data histories ufb and y using the ocid-identified system 
matrices A, B, C, D, identified observer gain G, and identified existing controller gain F. 
The reconstructed u^k) , y(k) stored as column matrices in ufb rec , y_rec, respectively, 
are computed from the following equations. 

x (k + 1 ) = (A + GC)x(k) + (B + GD)[ Ufl ,(k) + «,(*)] - Gy(k) 

y(k) = Cx(k) + D[u /b (k) + u t (k)] 

u fb (k) = -Fx(k) 

Note that the above state equation equation is the same as the usual observer equation 
expressed in terms of the prediction error y(k) - y(k) 

x (k + 1) = Ax{k) + B^pik) + M t (*)] - G[y(k ) - y(*)] 

= Ax(k) + + M t (*)] - Gy(k) + G{Cx(k) + D^k) + m,(*)]} 

= (A + GC)x(k) + (B + GD)[ Ujb (k) + u t (k )] - Gy(k) 


Algorithm: 

The function is programmed using a Matlab function dlsim. 
See also: 
odd 
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Purpose: 


y_esti 


Compute estimated outputs using an identified observer. 

Synopsis: 

[t,yhat]=y_esu(a,b,c,d,G,u,y,dt,nptsJlag) 

Description: 

Any typical observer has the following form 

x(k + 1)= Ax(k ) + Bu(k) - G[y(k) - y(*)] 
y(k) = Cx{k) + Du(k) 

where x(k) is the estimate of the state x(k) and y(k) is the estimate of the output y(k). The 
system matrices A, B, C, D, and the gain matrix G may be identified from input data u(k) 
and output data y(k) using the function okid. The function 
[t,yhat]=y_esti(A,B,C,D,G,u,y,dt,nptsflag) returns a column vector, t, for the time period 
corresponding to the desired number of sample points ( npts ), and an npts x m matrix, yhat, 
for the estimated outputs. The flag is set to 1 for plotting. The user will be prompted with 
the desired number of sample points to be reconstructed. Plots will be given to show the 
comparison between the real test outputs and the estimated outputs. 

See also: 

okid, okid_b, okid_fb, y_pred 
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ypred 

Purpose: 

Compute predicted outputs using an identified system model. 

Synopsis: 

[t,yhat]=y_pTed(a,b,c,d,u,y,dt,nptsflag) 

Description: 

The identified discrete system with the sampling time, dt, has the following form 

x(k + \) = Ax(k) + Bu(k) 
y(k) = Cx(k ) + Du(k) 

where x(k) is the state vector. The system matrices A, B, C, D, and the gain matrix may be 
identified from input data u(k) and output data y(k) using the function okid. The function 
[t,yhat]=y _prcd(A,B ,C ,D ,u,y, dt, npts flag) returns a column vector, t, for the time period 
corresponding to the desired number of sample points ( npts ), and an npts x m matrix, yhat, 
for the reconstructed outputs. The outputs y(k) are used here for comparison with the 
reconstructed outputs. The flag is set to 1 for plotting. The user will be prompted with the 
desired number of sample points to be reconstructed. Plots will be given to show the 
comparison between the real test outputs and the reconstructed outputs. 


See also: 

okid, okid_b, okid_fb, y_esti 


79 



yucovar, yucovfb, yucovb 

Purpose: 

Compute the left and right correlation matrices for least squares identification problem. 
Synopsis: 

\ybarf, vbarf \ =yuco var(«,y,/? flag ) 

[ybarb,vbarb]=yucov_b(u,y,p) 

[ybarfybarf, ybarb, vbarb]=y\icovfb{u,y,p) 

Description: 

Given / sample points of a system with m outputs, r inputs, and an assumed system order 
of p*m, the input matrix u is of dimension / x r whereas the output matrix y is / x m. The 
flag is set to 1 for long histories. Let y fc be partitioned as 

y,=[y<0) y(l) y( 2) - y(/-l)]; 

where each measurement y(i)(i=l ,2,...) has the dimension mxl. The least squares 
problem for the forward observer Markov parameters is posed as 

v V T = Y V V T 
Iff 777 

The matrices y^V/ and V f V f T have the same structure as the cross and auto correlation as 
seen in the following 


y f v / = 


£ y(i)u T (i ) J)y(i + l)v T (i) — ^y(i + p)v T (i) 

.1=0 i = 0 i =0 


v,v/. 


I - 1 x 

I «0> 7 (l) 
1=0 

1-2 T 
1 v{i)u l (i+1) 

;=o 

I v(i)u l 0+2) 
i=0 


1-2 x 
X u(i+l)v‘ 0) 
i = 0 

1-2 T 
X v0)v / 0) 

i=0 

/-3 x 
I v0)v 7 0+1) 
<=0 


1-3 ... 

X u(i+2)v (i) 
1=0 

/-3 T 

X v0+lK 0) 

i=0 

1-3 T 
X v0)v 7 0) 
i=0 


/-p-1 


'T' 

i=0 


v(i)u T (i+p) £ v(i)v T (i+p- 1) £ v0> r 0+p-2) 

i=0 /=0 


i=0 

l-p-l 


u(i+p)v T (i) 


v(i+p-l)v T (i) 

i=0 

I v0+p-2)v ; 0) 

i '=0 


'T 1 

/=0 


v0)v 7 0) 


where 


v(0 = 


'«(/)' 

yd) 


; / = 0,1,...,/— 1 


The summations are performed using inner products to reduce computational time. 
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[ybarf, vbarf\=yucovar(u,y,p flag) returns yV f T and v f V/ in ybarf and vbarf 
respectively. 

For the backward approach, let y fc be partitioned as 

y b =[y( 0) y ( D y( 2) - y(t-p)], 

where each measuremen ty(i)(7 = 7 ,2 , . . .) has the dimension mxl. The least squares 
problem for the backward observer Markov parameters is posed as 

v V T = Y VV T 

y_ b v b , b v b v b 

The matrices yVj and V b V b T have the same structure as the cross and auto correlation as 
seen in the following 


yyj = 


Lb 


^y(i> T (i) ]£y (i)v T (i + p) — £y(i)v r (i + 1) 
1=0 1=0 i=0 


I u(i)u T d) 
<=0 

l—p T 

I v(i+p)u r d) 


v(i+p-\)u T (i) 

1=0 

V v(i+l)« T (i) 

. i=0 



'l u(i)v T (i+p ) if «(i')v r (i+p-l) 

^ i=0 . i=0 

X v(i+p)v 7 (i+p) Z v(/+p)v r (i+/»-l) 

{ i= 0 i i= 0 

iT v(i+ p-\)v ^ (/+ p) Z v(i+p-l)v^ (i+p-l) 
i= 0 . i=0 

jf v(i+l)v^(i+p-1) X v(i+l)v^(i+p-l) 
i=0 i=0 


V u(i)v r (i+l) 
i=0 

/ *i+p)v r (i+l) 
i=0 

I - P T 

Z v{i+p-\)v l (i+l) 

i=0 

l~P T 

I v(i+l)v y (i+1) 
i=0 


The summations are performed using inner products to reduce computational time. 

[ybarb, vfrar7>]=yucov_b(M,y,p) returns y^V/ and V^V k T in ybarb and vbarb respectively. 
Note that there is no flag in this function! 

The strong similarity among these matrices y f V f T , y^l"/, V^/ f T and KK T suggests the 
possibility of simultaneously computing them in one single function. 

[ybarf, vbarf, ybarb, vbarb]=yucovfb(u,y,p) returns y f V f T , V f V f T > y^J an d V b V b T in 
ybarf, vbarf, ybarb, and vbarb, respectively. 



Example: 


y=[5 6 7 8 9]'; 
p=2,flag=0; 

tybar,ubar]=yucovar(u,y,p,flag) 
ybar == 

80 50 200 26 146 

ubar = 

30 20 70 11 56 

20 14 44 8 38 

70 44 174 23 128 

11 8 23 5 20 

56 38 128 20 110 

See also: 

arx_bat, arx_b, arx_fb 
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yycovar 


Purpose: 

Compute the left and right output residual correlation matrices. 

% 

Synopsis: 

\yvt,vvt\=yycovax(y,p,iexp) 

Description: 

Given / sample points of a system with m output residuals (the stochastic part of the 
system) and an assumed system order of p*m, the output residuals y is / x m. The flag is 
set to 1 for long histories. Let y be partitioned as 

y = [y(0) y(l) y( 2) - y(l- 1)]; 

where each y(i)(i=I ,2,...) has the dimension mxl. The least squares problem for the 
state estimator Markov parameters which whiten the residual is posed as 

yV T = YW t 

the matrices yV T and W T have the same structure as the cross and auto correlation as seen 
in the following 


yV T = 


£y(/ + l)/0‘) Xy0' + 2)/(i) — Jy(i + p)/(0 


i = 0 


(=0 


W T = 


1-2 


2y(0/(0 

i= 0 

£y(/)/0+l) 


/-3 


£y(i +i)/(0 

i*0 

Xy(i)/(0 


1=0 

(-4 


i= o 
1—4 


5>(0/(* + 2) £y(iV(i + l) 


i=0 

l-p-2 


/= 0 
/-P-2 


i=0 


f-4 


I>0+2)/(0 

i=0 

X>(i + !)/(-) 

1=0 


i=0 


&0/0 + P) &0/0'+P-1) ^W)/0 + p-2) 
1=0 1=0 1=0 


i£y0' + p-i)/(0 

i=° 2 

l>0’ + P-2)y r (0 

i=0 

&i'+P-3)/(0 


i=0 


where 

y = [y(i) y(2) — y(/-i)] 

The summations are performed using inner products to reduce computational time. 
\yvt,vvi]=yycovar(y,p,iexp) returns yV T and W T respectively in yvt and vvr. 


83 


Example: Some random numbers are created to verify the computation of this function. 

randCnormal'); 

y=rand(2,l 1); 
v=fy(:,1:10)]; 
yl=y(:,2:ll); 
yl*v' 
ans = 

-4.6867e-01 2.4383e+00 
-1.0534e+00 -3.8444e+00 
v*v’ 
ans = 

2.4726e+00 -2.0860e+00 “ 

-2.0860e+00 1.3314e+01 

y=y'; 

[ ybar,ubar]=yycovar(y, 1 ,0) 
ybar = 

-4.6867e-01 2.4383e+00 
-1.0534e+00 -3.8444e+00 
ubar = 

2.4726e+00 -2.0860e+00 
-2.0860e+00 1.3314e+01 

See also: 

k_abcd, arx_bat 
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